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I .   INTRODUCTION 

Biological  processes  have  been  operated  as  continuous  flow 
systems  for  many  years.   However,  development  of  mathematical 
models  for  such  systems  has  been  more  recent  (since  1950)  with 
significant  early  contributions  by  Novick  and  Szilard  (1), 
Monod  (2),  and  Herbert  (3).   Most  of  the  earlier  work  was  devoted 
to  the  analysis  of  pure  bacterial  cultures  in  single-stage  reactors. 

More  recently,  there  has  been  much  interest  in  the  analysis 
of  biological  waste  water  treatment  (^,5i6).   These  works  again 
assumed  that  the  process  Involved  a  single  culture. 

However,  anaerobic  digestion  is  a  complex  process  and  in 
the  recent  past,  growing  attention  has  been  given  to  it  from 
both  the  theoretical  and  experimental  points  of  view.   Malina  and 
HcCarty  (7,10)  refer  to  anaerobic  digestion  as  a  complex  two- 
step  process  involving  various  intermediate  chemical  species  and 
several  types  of  organisms.   This  mechanism  is  now  widely  used. 
Willlmon  and  Andrews  (9)  carried  out  experimental  work  using  a 
single  stream  system  under  various  operating  conditions.  These 
authors  (8)  have  given  a  mathematical  formulation  of  the  kinetic 
model  of  the  anaerobic  process  which  allows  them  to  simulate 
one-stage  and  two  stage  processes.   Finally,  Pfeffer  (13)  has 
emphasized  the  advantages  of  a  contact  process  (i.  e.  including 
recycle  of  organisms)  as  contrasted  with  a  conventional  system. 

In  this  work  a  mathematical  formulation  employing  the 
kinetic  model  of  Willlmon  and  Andrews  Is  used  to  simulate 
conventions?.!  and   contact  anaerobic  processes  consisting  of 


tt/o  stages.   In  addition,  by  adjoining  an  economic  model  to  the 
process  model,  the  process  is  optimized  for  various  values  of  the 
recycle  ratio. 

The  results  obtained  by  this  approach  must  be  balanced  with 
engineering  judgment  and  experience.   The  kinetic  and  economic 
models  can  only  approximate  reality.   Nevertheless,  this  study 
may  yield  a  better  understanding  of  the  process  and  a  more 
efficient  industrial  application. 

11  •   BASIC  .HgLATIOKSillPS 

1.   Mechanism  of  anaerobic  fermentation  (7) 
The  anaerobic  treatment  of  waste  water  involves  several 
microbial  species  which  carry  out  numerous  biochemical  and 
microbiological  reactions  (see  Figure  1).   The  process  yields 
carbon  dioxide  C02,  methane  CHjj,,  and  reduced  organic  molecules 
(HpS,  etc...).   The  bacterial  population  consists  of  facultative 
organisms  which  tolerate  small  amounts  of  dissolved  oxygen  and 
anaerobic,  bacteria. 

The  anaerobic  fermentation  process  is  a  sequential  one 
including  two  distinct  steps,  which  are  the  "acid  fermentation" 
step  and  the  "methane  fermentation"  step.   During  the  "acid 
fermentation"  step,  "acid  producing  bacteria"  break  complex  organic 
Compounds  down  to  simpler  organic  structures,  as  bacterial  growth 
takes  place.   The  principal  intermediate  compounds  resulting  from 
"acid  fermentation"  are  volatile  acids,  i.e.  short-chain  carboxylic 
acids  (C,  to  C;;).   These  volatile  acids  provide  substrate  for 
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the  "methane-froming"  bacteria.   During  "methane  fermentation", 
the  organic  acids  produced  during  the  "acid  fermentation"  step 
ere  converted  into  carbon  dioxid_e  and  methane.   These  bacteria 
are  substrate  specific,  i.  e.  each  of  these  ferments  only  a 
small  group  of  intermediate  compounds.  This  fact  has  also  been 
recognized  by  Hillimon  and  Andrews  (9)  in  their  experimental 
work.   Thus  the  stabilization  of  all  intermediates  necessarily 
Involves  several  cultures.   Figure  2  shows  the  significance  of 
acetic  and  propionic  acids  as  intermediate  products  (10). 

We  shall  implement  this  scheme  in  the  kinetic  model  by  using 
the  following  monenclature. 

S  =  raw  material  to  be  converted, 

A  =  acid  producing  bacteria, 

R  =  intermediate  product.   It  undergoes  fast  conversion  by 
methane-producing  bacteria, 

U  =  intermediate  product.   It  undergoes  slow  conversion  by 
methane-producing  bacteria, 

B  m   methane  producing  bacteria  fermenting  H, 

C  m   methane  producing  bacteria  fermenting  U, 

P  ■■  final  product. 

Figure  3  is  a  schematic  representation  of  the  mixed  culture 
model  which  is  assumed  in  this  study.   Although  there  are  a  large 
number  of  intermediates  and  microbial  species,  we  shall  assume 
that  the  system  can  be  adequately  represented  by  two  intermediates, 
R  end  *J,  and  three  microbial  species,  A,  B  and  C.   Figure  k   shows 
the  system  of  two  completely  mixed  tanks  In  which  the  porcess  is 
carried  out.   Such  a  system,  when  it  does  not  include  a  recycle 
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Fig.  2     Pathways   in   methane  fermentation   of 
complex    wastes   (10). 
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stream  of  organisms,  is  referred  to  as  a  "conventional  system" 
(10).  When  it  includes  a  recycle  stream  of  organisms,  it  is  a 
"contact  process"  (Figure  5)» 

2.   Kinetic  model 

Several  assumptions  must  be  made  in  developing  the  kinetic 
model  for  this  mixed  culture  system. 

1.  Environmental  conditions  are  such  that  acid  fermentation 
occurs  only  in  digestor  1  and  methane  fermentation  only  in 
digester  2. 

2.  Both  digesters  are  completely  mixed. 

3.  The  effect  of  endogeneous  respiration  or  organism 
deday  can  he  neglected. 

h.     Product  fermentation  is  directly  related  to  growth. 

5.  Isothermal  conditions  are  assumed  and  thus  there  is  no 
temperature  effect. 

6.  Monod's  function  describes  the  growth  rate,  1.  e. 

kSX 
rx  ~  Ks  +  X 

where 

r  =  rate  of  production  of  organism, 

k  -  maximum  specif io  growth  rate, 

K„  =  saturation  constant, 
3 

S  =  substrate  concentration, 
X  c  organism  concentration. 

For  the  sake  of  clearness,  we  shall  denote  the  various 
convert-ions  involved  in  the  kinetic  scheme  aa  follows  t 
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A 
(Al)   S  -  R 

A 
(A2)   S  -  U 

B 

(B)  R  -  P 

C 

(C)  U  -  P 

Consider  conversions  (Al)  and  (A2) ,  where  S  is  metabolized 

by  organism  A  yielding  products  R  and  U  and  additional  cells  A. 

The  rates  of  utilization  and  creation  of  the  various  species 

involved  are  directly  ralated  by  the  yield  factors.   The  same 

holds  for  conversions  (B)  and  (C).   For  instance,  yield  factors 

Y,  ,„  and  Yr,/q  are  defined  as 
A/3      va 

y       rate  of  formation  of  A 
A/3  ""  rate  of1  utilization  of  3 

v     rate  of  formation  of  B  by  reaotion  (Al) 

R/S  =        rate  "of  utilization  of  ~3 

The  various  yield  factors  are  shown  in  Figure  6. 

la  the  first  stage,  substrate  S  is  consumed  by  organism  A, 
giTing  organic  acids  R  and  U  as  products.  Using  Konod's  model 
for  growth  of  A  gives 

rA  =  J^  (1) 

From  this  expression  and  the  yield  factor  X»/g,  we  obtain  for  the 
consumption  of  substrate  (-rt, ) 
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Since  the  rate  of  production  of  organic  acids  H  and  U  is  directl.v 
related  to  growth 


\  -  -  YH/3  rs  -  ijj  *A  (3) 

*Vl  -  -  Yu/s  rs  -  r^i  rA  c*) 

In  the  second  stage,  intermediate  R  is  consumed  by  organism 

B  and  intermediate  U  is  consumed  by  organism  C  with  product  P 

(methane)  being  produced.   The  kinetic  models  for  growth  of 

organism  E  and  consumption  of  R  are,  respectively  (the  consumption 

of  R  is  -  r   ) 
R2 

r  _  W_  (5) 

XB  "  KR  +  R2 

kBR2B  (6) 

rR2  "  "  V?KH  +  K2) 

Likewise,  for  organism  C  and  intermediate  U,  we  can  write 

-     _W_  (?) 

0   -   Ky   +   U, 


*tl   = 


U2  =  "  *C/UiKU  +  °"2J 


(8) 


Since  product  P  is  produced  by  the  fermentation  of  both  B  and  C, 
the  expression  for  product  yield  is  of  the  form 

v     _  „  v  .  v   _  v    v  ( P ) 

rp  =  -  rp/B  r^   .p/u  r^  v.  ; 

or  in  terms  of  organisms  B  and  C 
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r  =  p&  rB  ♦  &2  r  do) 

p     yb/b   B     sc/u   c 

In  the  next  section,  these  equations  will  be  combined  With 
the  material  balances  at  each  stage  for  each  species  involved  in 
the  process.   This  will  yield  the  performance  equations  of  the 
system. 

3.   Performance  equations 

The  material  balance  equations  are  derived  first  for  a 
conventional  system  (see  Figure  k ) .   Then  they  will  be  developed 
for  an  anaerobic  contact  process  (see  Figure:  5) 

A.     Conventional  system. 

At  the  first  stage,  the  four  species  involved  are  organism 
A  and  organlcs  S,  H  and  U.  A  material  balance  for  organism  A 
yields  ' 

-A  ♦  rA  6;L  -  0 


A  =  r.  9.  (11) 

A  1 

Substituting  r  from  equation  (1),  and  solving  for  S.  ,  \ie   have 

si  -  e,kA  -  i  UV-J 

A  substrate  S  material  balance  gives 


Vi 


_rs8i  (13) 
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Dividing  equation  (11)  by  equation  (13)  gives 

A  .  (8   -  S  )  (  -  I*) 

S 
r 

But,  by  definition  -  ~A  =  YA/r,.   Therefore  we  have 
rg    A/3 

A  .  '^(S,  -  Sx)  (14) 

A  substrate  R  material  balance  for  the  first  stage  yields 

-  Ri  ♦  V1  =  ° 


(15) 
Dividing  equation  (15)  by  equation  (13)  gives 
h   «  (80  -  SlJ  (-  A) 


By  definition 

rP-l 

rg     -  XR/S 

Tills  gives 

Rl  n   YR/S(S0    " 

■  s. 

X 

(16) 


In  the  present  kinetic  scheme,  the  processes  of  formation  of  R 
and.  U  are  formally  identical.   Therefore,  we  nan  write 

"a  -  -u/s(so  -  3i>  <1?) 

In  the  second  stage,  five  species  are  involved,  namely  B,  C, 
R,  U  and  P.   Vie  can  write  the  steady  state  material  balances  for 
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each  component  as  follows!   For  organism  B  we   obtain 

-  B  +  rBe2  =  0  (18) 

or 

E -2 09) 

K2  -  kEe  -  1  v  y' 

Consider,  the  consumption  of  substrate  U  by  organism  C. 
The  kinetic  relationship  governing  this  growth  process  is  formally 
the  same  as  for  the  reaction  of  organism  B  on  substrate  R  [Equa- 
tions (5)  and  (7)].   Moreover,  the  material  balance  for  organism 
C  around  the  second  stage  has  the  sace  form  as  equation  (IS), 
As  a  result,  we  can  write 

The  material  balance  at  the  second  stage  for  intermediate 
R  is 

El  "  h   +  rR292  -  °  (21) 

or 

B,  -  B2  =  -  rR  e2  (22) 

Combining  this  with  equation  (18),  we  get 

B  .  Si.   (rx  -  r2)  (23) 


2 
Taking  account  of  equations  (5)  and  (6),  we  have 
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Due  to  the  similarity  of  the  reaction  of  B  on  R,  and  of  the 
reaction  of  C  on  U,  we  can  write  immediately 

c  ■  Yc/u(TJi  -  V  <?-5> 

The  steady  state  material  balance  for  product  P  is 

-  P  +**£e2  =  0  (26) 

Substituting  the  right  hand  side  of  equation  (10)  for  r  we  obtain 

P  „  _£Z?  r  o     +  -£Z2  re, 
YB/R  b  2   ~C/V     C  2 

Taking  account  of  equation  (IS)  and  the  corresponding  relationship 

C  ~  r  S  we  obtain  after  substitution 
C  2 


P/R    D  .  XP/U 


-  +  ^  .  C   ,  (2?) 

*B/R       XC/U 

Equation  (12),  (1^),  (16),  and  (17)  describe  the.  performance 
of  the  first  stage  while  equations  (19),  (20),  (2'0i  (25),  and 
(2?)  describe  the  behavior  in  the  second  stage  for  a  two  staye 
conventional  process.   The  corresponding  set  of  equations  for  a 
contact  process  which  includes  a  recycle  stream  of  organisms  will 
be  derived  next, 

B,   System  with  recycle. 

Since  the  model  for  growth  used  in  this  study  assumes  that 
organism  A  grows  only  In  stage  one  and  organisms  B  and  C  grow 
only  In  stage  two,  any  attempts  to  concentrate  the  organisms  by 
means  of  recycle  must  be  carried  out  by  feeding  the  organisms 
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leaving  a  stage  back  to  that  same  stage.   The  methane  producing 
bacteria  represented  in  this  study  by  organisms  B  and  C  often 
have  a  smaller  maximum  specific  growth  rate  than  the  acid 
producing  organisms  (organism  A  in  this  study).   Because  of  this, 
recycle  of  organisms  B  and  C  will  be  considered  in  this  study  as 
shown  in  Figure  5>   Organisms  B  and  C  which  metabolize  the 
volatile  acids  in  the  second  stage,  are  settled  in  a  clarifier  at 
the  outlet  of  this  stage  (Figure  5).   A  fraction  r  of  the  sludge 
leaving  the  clarifier  is  recycled  to  the  second  stage.   Such  a 
system  is  referred  to  as  an  "anaerobic  contact  process"  (10).   It 
is  assumed  that  there  is  no  endogoueous  decay  of  organisms. 
All  streams  after  the  second  stage  are  supposed  to  be  equally 
concentrated  in  oiganlos  S,  H  and  U.   The  recycle  stream  does  not 
contain  product  P.   The  treated  effluent  does  not  contain  signifi- 
cant quantities  of  organisms  B  or  C  (Figure  5). 
Let 

q  =  fraction  of  the  sludge  sent  to  waste  disposal  section. 

p  =  clarifier  efficiency. 

r  =  recycle  ratio. 

The  kinitic  scheme  given  in  Section  2  above  remains  valid 
for  the  present  system,  i.  e.  the  rates  of  reaction  of  che  various 
specie»  involved  in  the  process  are  given  by  equations  (1)  through 
(10).   Similarly,  for  the  first  stage,  there  is  no  modification 
to  equations  (12),  (1*0,  (16),  and  (1?)  needed,  because  there  is 
no  change  in  performance  at  this  stage. 

The  performance  equations  for  the  second  stage  are  obtained 
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from  the  steady  state  material  balances  for  species  B,  C,  H,  U, 
end  P. 

The  amount  of  organism  B  entering  the  second  stage  is  now 
(qr)  ( p B ) .   Moreover,  the  flow  rate  through  this  stage  is 
q(l  +  r) .   Eence,  the  material  balance  for  B  can  be  written  as 

ft' 

qr?B  +  rBV2  -  lf(l  +  r)B  =  0  (28) 

where  rfi  is  given  by  equation  (5).   Substituting  this  expression 
for  rE  into  equation  (23)  and  dividing  by  q  yields 

r3B  +  P~2~   .  ST-  (1  +  r)B  E  0 


B 

"  Il2 

Re 

arr 

anging 

this 

gives 

kBR2 

h 

KR  + 

x1R2 

'   q 

1  +  (1  -,g)r  (29) 


By  definition  V"2  •     =  9   is  the  hydraulic  residence  time  at  the 
second  stage  and 

t 

g=  1  +  (1  -  §)r  ■  (30) 

is  a  dimensionless  factor  which  depends  on  g  and  r.   Then, 
equation  (ii.)  can  be  written  as 

kr  r^  =  * 

Solving  this  equation  for  H2  gives 

■a-  kB«2  -5: 
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Dividing  the  numerator  and  denominator  of  this  fraction  by 
yields 

where 

kB'  -  Vs  (32) 

Equation  (31)  has  the  same  form  as  equation  (19) .  where  k_  is 
replaced  by  k„'.  Note  that  equation  (28)  may  be  written,  after 
rearranging,  as 

rBe2=$B  (33) 

In  the  process  of  metabolization  of  specie  0  by  organisms 
C,  the  underlying  kinetic  scheme  is  the  same  as  for  the  reaction 
of  R  and  3.   Therefore,  by  analogy,  we  obtain 

where 

V  -  kc/^  (35) 

The  amount  of  R  carried  by  the  recycle  stream  is  (rq)R0, 
The  amount  of  R  leaving  the  second  stage  is  q(l  +  r)R?.   Then, 
the  steady  state  material  balance  of  this  compound  is 

I'qK,  +  qRl  +  rp  V,  -  q(l  +  r)R,  a  0  (36) 

Searrangins  equation  (36) ,  we  obtain 
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/ 

rR262  "  R2  "  h                                                                                                                <37) 

Dividing  equation  (33)  by  equation  (37)  yields 

3'b      b  1 

rH2  =  \   ~   Rl 

Prom  equations  (5)  and  (6) 

rH2  "   B/R 

Substituting  this  value  into  the  precseding  equation  yields 

B .  lm  (h  -  b j 

^     1    2 

or 

B  =  *3/H(Hl  "  H2)        ^                          (33) 

where 

1  B/R  -   |  (3; 

Equation  (33)  has  the  same  form  as  equation  (23)  where  ¥._,/,,  is 
replaced  by  Y'B/R- 

The  steady  state  material  balance  for  organic S  U  can  be 
obtained  by  analogy  to  the  case  of  organics  E.  Hence  we  can 
write  Immediately 


'••  -  *W°i  -  V  {h0) 


wnere 


*'c/u  -  -f-  <*" 
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For  product  P,  the  steady  stats  material  balance  Is 


P  (1  +  r)q  t  rp  V2  .  0 


Substituting  rp  as  given  by   equation   (10),   we  obtain 

P        VhV2      +    VcS  (^j 

According  to  equation  (33)  v'o  have 

rBe2  „$B 

By  analogy  we  can  write 

rO  =  f  C 
C  2   * 

Substituting  these  expressions  into  equation  {k})    yields 

*  fcir?r5  YB/R  +  (l  ♦-?)  yc/u  «     ' 

Let 

-  P/B   1  +  r 


1  P/U  "  1  -;-  r 


C*5) 


Then,  taking  account  of  equations  (39),  ('*!),  arid  (45),  equation 
(n-'i  i  can  be  written  as 
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1  3/R       X  C/U 

Performance  equations  (3D.  (3*)i  (30),  (40)  and  (46)  have  the 
same  form  as  the  corresponding  equations  for  the  conventional 
process,  i.  e.  equations  (10),  (20),  (24),  (25)  and  (27) 
respectively.   Therefore,  provided  the  constants  kg,  k  ,  Y  .  , 
X „/,,,  Yp/R  and  Yp/n  are  transformed  according  tc  equations  (32), 
(35),  (39),  ('H)  and  (45)  respectively,  the  simulation  of  the 
process  can  be  carried,  out  using  the  same  computer  progratu  as 
for  the  conventional  process. 

An  important  parameter  for  a  contact  process  is  the  solid 
retention  time.   According  to  KcCarty  (10),  the  solid  retention 
time  (SHT)  is  defined  as 

S R'T   _§.H-S o ended  solids  in  the  system _ 
**  Suspended  "solids  removed  per  day 

We  shall  apply  this  definition  to  the  organisms  E  in  the  second 
stage.   The  mass  of  organisms  B  in  this  stage  is  V_B.   The 
quantity  of  B  removed  per  day  is  (  q)(3B).   Then 

V2B 
SlUB  "  USJb 


V?    1 
OK1B  -  q   *  wg  ^'> 


where  SRT„  is  the  solid  retention  time  of  organisms  B  in  the 
second  stage.   Furthermore  the  material  balance  for  organisms 
around  tfte  clarifyer  can  be  expressed  as 
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* 

(1  +  r)  qB  =  (r  +  n)    qpB 

which  yislds,  after  simplifying  and  rearranging, 

pw.  1  +  (1  -  p.)  r  (48) 

Then,  according  to  equations  (30),  (47) ,  and  C+8),  SRT  can  be 

o 

expressed  as 

SRTD  =  S_ 
B   % 

Note  that  equation  (^8)  implies 

i  +  (1  -  g)  r  >  0 
or 


vihere  ' 

rmax  =  g  -  1 

When  r  =  r   ,  equation  (^3)  shows  that 
max 

g~>=  0 
If  p  ^  0,  then  oo  =  0,  i.  e.  all  available  organisms  B  a-  recycled. 

C.   Concept  of  wash-out  time. 

Since  microbial  growth  is  an  autocatalytic  process, 
organisms  must  be  present  for  growth  to  occur.  If  the  flow  rate 
is  too  large  relative  to  the  reproduction  rate  the  organisms 
may  not  be  able  to  reproduce  fast  enough  and  wash-out  of  organisms 
may  result.   It  is  important  to  determine  which  residence  times 
will  continuously  support  a  growing  culture  and  which  ones  will 
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result  in  wash  out  and  no  growth.  At  first  we  shall  consider 
the  first  stage.   The  function  which  describes  the  variations  of 
outlet  concentration  of  substrate  S  is  hyperbolic  as  given  by 
equation  (12). 

KS 
Sl  "  *Jk/.  1  <12> 

There  is  a  mathematical  limit  for  9,  for  which  S-,  -  ro.   This 

1/k..  Obviously,  it  has  no  physical  signlfi- 

for  which 
equation  (12)  has  a  physical  significance  is  where 

si=so 

This  is  to  say  that,  when  9,  =  C  -,,  wash-out  occurs  for  ojg&nism 

A  and  there  is  no  conversion  of  substrate  S.   Equations  (16), 

(1?),  and  (14-)  show  that  under  this  condition  R,  =  U.  =  A  =  0. 

These  results  are  illustrated  in  Figure  ?.   When  6,  >  9  , . 

1     Wl 

Organism  A  can  grow  and  S  is  converted  to  produce  R  and  U. 
Equations  (12),  (1*0  ,  (16),  (17)  show  that,  under  this  condition, 
we  have 

fA  >  0 
H  >  0 
°1>  ° 

At  the  inlet,  of  the  second  stage  the  concentration  of  S,  A, 
B,  and  U  are  respectively  S^  ,  A,  f^,  and  0_.   Consider  equation 
(19)  which  gives  H  ,  the  outlet  concentration  of  H,  as  a 
function  of  residence  time  B      in  the  second  stage.   As  in  the 


?5 


s5d;s    jsjjj   ayj  oi  suoudj-jusojoo 
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previous  case  the  relative  size  of  the  values  of  the  maximum 

specific  growth  rate  k  and  the  dilution  rate  9  Will  determine 
B  2 

whether  growth  of  organism  B  can  occur  in  the  second  tank.   E  is 

converted  only  if  organism  B  is  present  in  which  case  the  inlet 

concentration  R,  exceeds  the  outlet  concentration  R  .   Examining 
1  2 

equation  (19)  then  shows  that  there  exists  a  residence  time  8. 

wz 

for  which  the  following  equation  holds. 

B2  -  Bj  »  0  {kQ) 

At  8   ,  equation  (23)  shows  that  B  =  0.   When  0  >  6   ,  R 
"2  2    w2 

begins  to  be  converted  and  B  to  grow,  giving  product  P. 

A  similar  reasoning  applies  to  U  and  C.   The  former  is 
converted  and  the  latter  grows  only  if  6  >  6   ,  where 


v3 


the  solution  of  equation  (20)  for  e   if 


u2  -  V1  (50) 

At  6  „,  equation  (25)  shows  that  C  =  0.   When  0„  >  8  „,  B  bagins 
w3  2    w3 

to  be  converted,  and  C  to  grow,  yielding  product  P. 

It  has  been  noted  (8)  that  in  the  second  stage,  wash-out 

usually  occurs  first  for  the  organism  which  has  the  smallest 

maximum  specific  growth  rate.   Since  k  <  k  ,  C  is  the  slowest 

C     B 

growing  organism.   Therefore,  the  following  inequality  is 
usually  satisfied. 

Ttms,  the  condition  which  must  be  satisfied  in  order  for  organism 
A  to  grow  in  the  first  stage  is 


2? 


61  *   \l  (52) 

Similarly,  the  condition  for  the  growth  of  organisms  B  and  C  in 
the  second  stage  is 

6  *  6  ->  '  (53) 

The  inequalities  given  by  equations  (52)  and  (53)  will  be  referred 
to  as  "conditions  of  feasibility". 

4.   Economic  model. 


We  shall  consider  a  cost  of  treatment  consisting  of  two 
terms.   The  first  term  is  the  total  hydraulic  residence  time 
6  +  6„  of  the  feed  stream  q  in  the  system.   This  term  represents 
approximately  the  cost  of  the  digesters  in  terms  of  volume. 
Therefore,  it  also  represents  roughly  the  fixed  costs  of  treatment 
for  a  given  plant  and  a  given  flow  rate.   In  addition,  the  degree 
of  removal  of  organics  from  the  feed  stream  has  an  important 
economic  significance.   In  fact,  the  stream  leaving  tile  clarlfier 
can  either  be  discharged  into  a  receiving  stream  cr  treated 
further  in  order  to  achieve  a  better  stabilisation  of  its  organic 
content.   In  both  cases,  we  may  incur  a  penalization  which  is  a 
function  of  the  degree  of  treatment.   In  this  work,  the  penaliza- 
tion term  takes  the  form  z„  S,  +  Zj,  VU   +  z,,  0„  penalization 
where  z  .  zQ,  and  z,,  are  constants.   This  term  must  also  be 

fa      Q  U 

expressed  in  terms  of  residence  time.   Thus,  the  dimension  of  the 
constants  z  ,  z  ,  and  z      is  time/concentration. 

III.   ANALYSIS  OF  SYSTEM 
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In  this  section,  we  shall  determine  the  wash-out  times, 

9   ,  S,  „ ,  S  _,  of  the  system  using  the  concept  of  Mash-oat  stated 
wl   *•  2   wj 

in  Section  II, 3 •  Then,  we  shall  prove  that  the  wash-out  steady- 
state  is  a  solution  for  the  steady  state  problem  provided  that 
the  feed  stream  is  free  of  organisms. 

1 .   Determination  of  wash-out  times. 
Organism  At 

3y  definition  of  6  , ,  we  have 
*  wl 

W  ■  so 

Because  of  equation  (12),  this  gives 

KS 

S, 


e  ,  k„  -  i   o 

wl   A 

and 

K   and  S„  are  positive.   Thus,  equation  (5k)    implies 

6   >  -I-  (55) 

wl   k , 
A 

Moreover,  when  S   n  5  ,  equation  (1A)  shows  that  A  «=  0.   When 
wash-out  occurs  in  the  first  stage,  the  concentration  of  organism 
A  is  zero. 


Organise  3i 

The  corresponding  wash-out   tslme,    say  0   „,   occurs  when 

w2 

W  -  Ri  (56) 
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Substituting  from  equation  (19) »  we  have 

a. 


kBew2-1"'1 


Then 


•  9w2  -  t(1  +  if  (57) 

Again,  we  notice  that  substituting  equation  (56)  into  equation 
(24)  leads  to  B  =  0. 

Organ! sq  C i 

6   corresponds  to  the  condition  that 
w3 

vy  -  ui  (58) 

Substituting  from  equation  (20)  we  have 
KU 


k„  e  ,  -  1  -  vi 

c  w3 


Solving  for  0  _  yields 
w3 


K 
tt(1  *  *«2)  (59) 


N   *3   V"  ■   ux 

This  equation  implies  that 


0  ,  >  -/-  (60) 


Equation  (5'0  shows  that  e  ,  does  not  depend  on  the  residence 
time  of  any  of  the  stages.   However,  the  expressions  for  0   and 
6  ,  contain  the  terms  H1  and  U,  respectively,  which  are  dependent 
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on  the  residence  time  of  the  first  stage.   As  a  result,  the 
wash-out  times  of  organisms  B  and  C  depend,  on  the  residence  time 
in  the  first  tank.   Figure  8  illustrates  the  concept  of  wash- 
out time  of  organism  A  in  the  first  stage.   The  dimenslonless 
concentrations  S./Sq,  A/Sq,  H^Sq,  and  U./Sq  as  obtained  from 
equations  (12),  (l^),  (16),  and  (17)  respectively  are  plotted 
versus  the  residence  time  0-,  ranging  from  0  to  2  days.   When 
8.  <  6wl.  si/so  =  1  and  A/so  -  °»  no  conversion  takes  place  and 
wash-out  of  organism  A  occurs.   When  S    >    3  ,»  the  concentration 
of  substrate  S  falls,  indicating  that  S  is  converted,  while  the 
concentration  of  organism  A  rises  due  to  its  growth. 

2 •   Steady  state  wash-out  time . 

Let  us  consider  the  first  stage.  If  we  assume  an  initial 
concentration  A  of  organisms  in 'the  feed  stream,  the  material 
balance  for  organism  A  can  be  written  as 

AQ  -  A  =  -  rA  Bj  (61) 

For  substrate  S,    it   Is 

3Q  -  sx  .  -ra  h  (7) 

Dividing   equation   (61)    by   equation   (7)    and  noting   that   ~r  /r     - 

A      3 

X .  ,„,   we   obtain 
A/S 

A0    -   A  "   "   YA/S(S0    "   V  <62) 


which   eau   be  rearranged    to  give 
km  A     -:•   I         <3      -   3    ) 
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Substituting  this  expression  into  equation  (1)  giving  r  yields 

A 

kAsl .  ..  r.       .    «        i*  a   \t  (63) 


rA  "  V^I  U°   +   WS0   "   Sl)] 


Substituting   equations    (62)    and    (63)    into   equation   (6l)   yields 

"   YA/S(S0    -    Sl)    =    -  J~-s~  CA0   +   YA/S(S0    -    Sl)]  (M) 

By  introducing  the  condition  AQ  =  0,  equation  (6^)  can  be 
written  as 

(so-si)(l--|A.Jl),o  (65) 

This  is  an  equation  of  second  order,  which  gives  the  steady  state 
concentration  S-.   Obviously,  one  of  the  roots  is  S  ^  S„,   This 
is  precisely  the  condition  for  wash-out  of  organism  A  from  the 
first  stage.   Another  root  is  given  by 

1   k*L_el_si  _ 
KS  +  31  ' 

Solving  for  3.  yields 

1   « 
which  is  identical  to  equation  (12).   Thus  this  approach  shows 
that,  in  addition  to  the  steady  state  described  by  the  performance 
equations  given  in  Section  II. 3.,  there  is  a  second  solution 
which  can  be  referred  to  as  "wash-out  steady  state"  for  the 
first  stage.   By  analogy  the  same  conclusions  can  be  reached  for 
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the  second  stage. 

The  analysis  given  above  yields  an  important  conclusion  as 
to  the  operation  of  the  system.   Indeed,  let  us  assume  that 
wash-out  of  organism  A  occured  in  the  first  stage  due  to  flooding. 
Then,  for  this  stage,  the  stable  steady  state  solution  is  given 
by  the  first  root  of  equation  (65),  i.e.  S  =  S„.   This  implies 
that  the  feed  stream  flow  through  the  first  stage  without  tinder- 
going  any  biochemical  conversion.   However,  if  a  residence  time 
6   is  established  such  that  6,  >  9k1  and  if  some  organism  A  is 
introduced  into  the  first  stage,  this  organism  can  grow  utilizing 
substrate  S.   After  a  transient  period,  the  system  reaches  a 
stable  steady  state  described  by  the  second  root  of  equation  (65) 
which  yields  equation  (12).   In  this  steady  state,  continuous 
organic  utilization  occurs,  which  is  the  condition  for  the  system 
to  operate  efficiently. 

IV.   OFTIHIZATIOK 

In  this  Section,  we  shall  determine  the  optimal  policy  for 
a  two-stage  continuous  anaerobic  digester  system  using  two 
different  approaches.   In  the  first  approach,  differential 
calculus  yields  an  analytical  expression  for  the  optimal  policy. 
In  the  second  approach,  the  analysis  of  the  problem  is  carried 
out  from  the  point  of  view  of  an  empirical  search  technique, 
namely  the  Simplex  technique  which  is  discussed  in  Appendix  IV. 

1 .  Formulation  of  the  problem 

According  to  the  economic  model  described  in  Section  II. 4, 
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the  objective  function  to  be  minimized  consists  of  two  parts. 

The  first  part  expresses  the  equipment  cost  (8  +  9_).   The  second 

pert  is  the  penalization  cost  z„   S»  +  z  R  +  z  U  for 

D   JL      R   2      U   2 

discharging  the  organics  remaining  in  the  effluent  stream.   Thus 
the  objective  function  J  can  be  expressed  as 

J  =  ex  +  e2  +  zg  $x  +  z    R2  +  zD  u2  (66) 

Substituting  equations  (12),  (19),  and  equation  (20)  into  the 
above  equation  gives 

Z„    K0  Zrv  K 


(6?) 


\   kA  "  l        kB  32  "  X    ■  kC  62  "  X 


Note  that  if  S  „  <  8~  <  6Tr0,  the  concentration  of  U  in  the 
w^  —   2    wj 

treated  effluent  is  U  .   Then  the  objective  function  WOttld  be 

z     K     '   z  K 
-•   -  -  t  z„  0, 


J  =  wi  +  H2  +  a.  k.  »  T  T  STe,  -  l  T  u  1 

1   A        B   2 

The  decision  variables  are  chosen  to  be  the  two  residence  times 

6  and  8,.   Since  there  is  no  equality  constraint,  the  optimlza- 
1      £- 

tion  problem  is  two-dimensional.   Equation  (6?)  shows  that  the 
optimization  problem  is  a  non-linear  one  with  the  decision 
variables  9   and  9   subject  to  the  conditions  of  feasibility. 

V  8wi  (52) 

2 •  Numerical  data . 

The  optimization  will  be  performed  with  the  following  values 
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of  physical  and  economical  parameters i 


so   - 

10.0  gm/1 

kA   = 

6.0  day"1 

Ks   ■ 

O.SO  gm/1 

YA/S  = 

0 .  25 

Y    - 
R/S  - 

O.ii-0 

Yu/s  = 

0.20 

kB   " 

0.50  day"1 

KR   = 

1.00  gm/1 

iB/R  = 

0.10 

Y 
P/R  - 

0.75 

kc   B 

0.25  day"1 

Ku   " 

1.50  em/i 

Yc/u  ** 

0.10 

Yp/u  = 

0.75 

zs   ■ 

1.00 

z    ■ 
R 

1.50 

z„   - 

1.50 

The  choice  of  these  values  for  zs,  z  ,  and  z  is  compatible  with 
the  two  conditions  of  feasibility  as  will  be  shown  in  Section 
IV.  3  below. 

3«   Analysis  by  differential  calculus. 

Lot  (L)  be  the  surface  described  by  equation  (67).   First, 
we  have  to  find  the  stationary  points  on  surface  (£).   The 
coordinates  of  such  points  satisfy  the  necessary  conditions 
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(68) 


and 

ft  "  °  (69> 

From  equation  (66)  we  have 

_?J    ,      ffl      dE2      dU? 

537  =  -1  +  zs  ae^^  +  zr  de1  *  zu  dex 

R  and  ti2  do  not  depend  on  6  .   Then  substituting  dS  /dS-,  yields 

zq  j:s  k„ 
1 §_  S   a   =  0  (?0) 

(kA  e1  -  l)2 


Solving  for  9..  yields 


In  order  for  6,  to  have  a  physical  significance,  it  has  to 
satisfy 

0,  >   0  .  >  X 
1    wl   ka 


Therefore,  the  expression  for  8.  is  necessarily 

6il  -  rA[1  +  ^PT^l  3  <72) 

Let  us  now  consider  equation  (69).   Prom  equation  (66) 
v  .  dR       dU 

i>8™  "  -  +  zr  de2  +  zu  de~ 
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dR  dU 

Substituting  into  this  equation  the  derivatives  — 0  and  — § 

a62  d02 
computed  from  equations  (19)  and  (20)  yields 

aj  _  ,     ZR  KR  kB       ZU  KU  kC 

>Q     ~   X  2  2  =  ° 

(kB  °2  '  1}  .    (kC  ?2  "  « 


+  — f2-_u_c   _  j  =  0  (?3) 


(kfl  e2  -  i)    (kc  e2  -  i) 


Let 


z  K„  k„       z„  K,.  k, 

f ( e2)  ,  — R_A_b_  +  _JL-L_£__  _.  i  (7-:) 

(kB  e2  "  1J     (kc  °2  "  1} 
Then,  equation  (7^)  becomes 

fi&2)   =  0.  (75) 

Equation  (73)  can  also  be  written  as 

ZRKRk3(kC    V    1)2+    ZUKUkC(kE   V    1] 

-   (kB  e2-  l)s(kc  9?_-  l)z  .  0 

This  equation  is  of  4th  order  and,  therefore,  equation  (73) 
has  four  roots. 

In  order  for  Gp  to  have  a  physical  significance,  it  has  to 
satisfy 

62  *   Sw3  >  i| 
It  can  ba  shorn  (see  Appendix  II)  that  the  largest  root  of 
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equation  (73)  f  say  6   ,  always  satisfies  this  condition.  It 

remains  to  prove  that  the  point  defined  by  the  coordinates  6 

1L 

and  9OT  is  a  minimum  on  response  surface  (2  ) ,  i.e.  that  the 
Hessian  matrix  of  the  objective  function  J  is  positive  definite 
at  this  point  (sea  Appendix  III). 

Equations  (72)  and  (75)  allovr  us  to  express  the  inequality 
constraints  (52)  and  (53)  in  terms  of  the  physical  and  economical 
parameters  of  the  system. 
The  first  Inequality  constraint  is 

6   >  9 
1L    wl 

Substituting  equations  (5't-)  and  (72)  in  this  equation  leads  to 
the  relationship 

i[i  +  vipc7^]>fci  +  ^i 

A  A       0 

which  must  be  satisfied  to  have  growth  in  the  first  tank. 
Simplification  of  the  above  relation  yields 


Jz     K  k  >  -2 
S   S   A   SQ 


Raising  this  to  the  second  power  and  rearranging  gives 

ZS  kA  (V3  "  h   >   °  (7S) 

Values  of  all  the  parameters  appearing  in  this  inequality  are 
positive.   It  foliov:s  from  inequality  given  by  equation  (76)  that 

i    > § 

kA  <S(/ 
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For  a  given  initial  substrate  concentration,  there  is  a 

lovrer  limit  for  penalization  z     in  order  for  the  first  stage  to 

s 

operate  above  the  wash-out  residence  time  and  under  optimal  con- 
ditions z„  -  1.0  satisfies  the  inequality  of  equation  (76). 
The  second  inequality  constraint  is 

e   >  e  .  (77) 

2L    »3 

where  9   is  the  largest  root  of  $f«2)  =  0  and  e^  is  given  by 
equation  (59).   Vie  know  that  0  ~  and  6    are  both  larger  than 
l/kr,  and  that  f(8?)  monotonously  decreases  vrhen  02  >  1/1.,,. 
Therefore,  the  inequality  given  by  equation  (77)  yields 

f(ev;3)  >?(e2L) 

But,  by  definition, 

<ne2Ij)  -  o 

Hence,  the  second  condition  for  feasibility  is 

f<ew3)  >  o 

Substituting  the  expression  for  f   from  equation  (7^)  yields 

_^H_1h13_^  +   ******      ,  -  1  >  0  (78) 

(kB  ew3  -  l)   (kc  ew3  -  i) 

This  inequality  is  linear  in  zp  and  e ,   In  addition,  6   depends 
on  6- .   It  can  thus  be  concluded  that  for  given  operating  con- 
ditions at  the  first  stage,  there  are  lower  limits  for  penaliza- 
tion coefficients  Z„  and  z„.   z_  =,  z     ^  1.50  satisfies  the 
It       u     H     U 
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Inequality  given  by  equation  (70). 

Finally,  the  coordinates  9   and  0   of  point  A  minimize  the 
objective  function  J  given  by  equation  (6?)  and  satisfy  the 
inequality  constraints  given  by  equations  (52)  and  (53).   Thus, 

8     -  0 
10PT  ~~  1L 

620PT  =  62L 

Point  A  corresponds  to  the  optimal  policy  on  response  surface  (£ ) . 
In  addition  to  point  A,  the  surface  (E'J  has  three  other  stationary 
points,  namely,  B,  C,  and  D  (see  Appendix  III).   B  and  D  are 
saddle  points,  C  is  a  maximum  point. 

k .     Analysis  by  empirical  search . 

As  stated  previously  the  objective  function  J  is  to  be 
minimized  under  the  two  inequality  constraints  8..  >  0  n  and 
8  >  6  _,  where  9  and  9  are  the  two  decision  variables. 

In  the  differential  calculus  approach,  this  constrained 
optimization  problem  is  solved  in  two  steps.   In  the  first  one, 
we  determine  the  stationary  points  of  the  response  surface  (  £ ) . 
In  the  second  one,  we  check  a  posteriori  that  the  stationary 
point  of  interest  (point  A)  satisfies  the  inequality  constraints 
given  above. 

The  empirical  search  technique,  Simplex,  is  particularly 
suitable  to  treat  at  once  such  a  constrained  optimization  problem 
(see  App  end  i  x  I V ) . 

By  Beans  of  these  two  techniques,  it  is  found  that  the 
optimal  residence  times  for  the  system  under  consideration  are 
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Ox  ■  0.455  days  and  62  =  7.183  days.   Figure  9  shows  a  plot  of 
the  dimensionless  concentrations  R_/S  ,  V„/SQ,    E/S  ,  C/3  and 
P/3  in  stage  tvjo  for  the  optimal  policy  for  stage  one.   Figure 
10  shows  the  plot  of  the  same  dimensionless  concentration  for  a 
contact  system  where  r  =  0.25' and  3  =  4.0.  These  figures 
illustrate  the  increased  stability  of  the  contact  process.   Indeed, 
the  wash-out  times  8   and  6  _  for  the  contact  process  are  0,629 
days  and  1.772  days  as  compared  to  2.515  days  and  7.037  days  for 
the  conventional  system.  Comparing  Figures  9  and  10  shows  also 
a  strong  increase  in  organism  concentrations  B  and  C  together 
with  a  more  efficient  reduction  of  the  organics  in  the  case  of 
the  contact  process. 

5.   Contours. 

When  J  is  given  a  fixed  value  J  equation  (67)  represents 
the  corresponding  contour.   Let  us  write  equation  (67)  In  the 
form 


where 


kB  o2  -  i   kc  e2  -  i 


(79) 


(80) 


Rearranging  equation  (79)  yields 

kA(01f  +  [(A  -  Jc)  kA  -  1]81  +  zg  Kg  +  Jo  -  A  .  0    (81) 

For  the  purpose  of  generation  of  the  contours,  let  y  be  the 
current  value  of  8„,  and  x..  and  x  the  current  values  of  the 
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two  roots  of  equation  (81).   Thus  the  procedure  to  get.  the 


Fix  y  (or  9  ). 

Compute  A  by  equation  (80). 

Substitute  A  and  solve  equation  (81)  for  x.  and  x  . 

The  two  points  (z.,  y)  and  (x  ,  y)  lie  on  the  contour  where 

the  obiective  function  has  the  value  J  .   By  allowing  y  to  vary 

c 

the  whole  contour  is  generated.  Thus  the  question  arises! 
what  is  the  permissible  range  for  y?  The  discriminant  Of 
quadratic  equation  (81)  is 

=  [(A  -  Jc)kA  -  if  -  4kA(zs  Ks  +  Jc  -  A) 

For  equation  (81)  to  have  two  real  roots,  has  to  be  positive. 
This  restricts  the  range  of  variation  of  A  and  9-  in  case  of  a 
maximum  or   minimum  for  response  surface  (E )  . 

The  contours  of  response  surface  (E )  around  the  stationary 
points  A,  B,  C  and  D  are  shown  on  Figures  11,  12,  13,  and  Ik, 

On  Figure  13,  the  contours  generated  by  this  procedure 
have  no  physical  significance  when  9w2  <  9?  <  9i  ,.   This  part  of 
the  contours  are  represented  by  the  dashed  lines.   In  deed,  when 
9   <  9  <  6  ,  wash-out  occurs  for  organism  C.   Thus,  the 
objective  function  is  given  by  the  following  equation  (see 
Section  IV.) 

J'  =  61  +  92  *  k/l^!  i  +  k/f ^  1  +  ZU  h  (82' 

where  U1  is  given  by 
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(12) 
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"   kA 

el   ■ 

-  1 

Equatl 

ons 

(12) 

and. 

(1?) 

yield 

ii 

V 

t< 

KS 

l  -  u/s  v  o   kA  e^^  -  I' 


Substituting  this  expression  of  U,  into  equation  (82)  gives 


z  K„ 

J-  =  9l  +  9  +  - 


B   2 


+  ZU  Vs(S0  -  i-4f'-T-T)  (83) 

When  J'  is  given  a  fixed  value  J'  ,  equation  (SO)  represents  the 
corresponding  contour.   Let  us  write  equation  (83)  in  the  form 


*'  =  e,  +  v^r- rr  ♦  zn  *„/«,  Sn  (85) 


where 

''Z   +  k^~9^~l  +  ZU  I0/fl  S0 
Rearranging  equation  (84)  yields 

kA(0l)-\[(A'  -  J'c)kA  -  l]Ga  +  H   Ks  -  Zu  Yu/3  Ks 

+  (J'c  -  A')  =  0  (86) 

The  contours  can  then  be  generated  using  the  procedure  previously 


50 


described.   They  are  shown  on  Figure  13  by  the  continuous 
lines  where  92  <  9   .   The  contours  of  the  physical  response 
surface  are  discontinuous  for  e_  =  6  „.   Indeed,  when  this 
equality  holds,  the  expression  of  the  objective  function  changes. 
This  explains  the  points  of  discontinuity  on  Figure  13. 

V.   RESULTS  AMD  DISCUSSION 

The  numerical  results  of  this  work  have  been  obtained  using 
the  program  MICU20  (see  Appendix  V).   The  main  design  variables 
for  the  conventional  process  and  a  contact  process  are  given 
below. 

For  the  conventional  process,  the  optimal  policy  is  S,  ■ 
0.455  days  and  8  =  7.183  days.   The  corresponding  concentra- 
tions are 
»  At  the  first  stagei 

S  =  0.289  gm/1 

A  =  2.4-28  gm/1 

B,  =  3.885  gm/1 

U1  =  1.91*2  gm/1 

-  At  trie  second  stagei 
Ro  =  O.386 
U2  =1.835 
3  =  0.350  gm/1 
C  =0.006  gm/1 
P  .  2.667  gm/1 

The  wash-out  times  are 
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6  ,  =  0.175  days 

V"  2-515deys.  . 

6w3  ■■  7.08°  days 

For  the  contact  process,  the  results  are  computed  for 
r  ■  0.25  and  B  ^   4.0.   The  optimal  policy  is  e  =  0 .^55  days 
and  9  =  2.64  days.   The  corresponding  concentrations  are 

-  At  the  first  stage 

Sx  =  0.289  gm/1 

A  =  2. 'I- 28  gm/1 

R  =  3.885  gm/1 

U  =  1.9*2  gm/1 

-  At  the  second  stage 

H  =,0.234  gm/1 

2 
U2  =.■  0.914  gm/1 

B  =  1.460  gm/1 

C  m   0.411  gm/1 

P  =  2.807  gm/1 

The  wash-out  times  are  • 

6wl  .  0.175  days 

6  „  -  0.629  days 
W2 

B  ,  -  1.772  days 
The  solid  retention  time  in  the  second  stage  is 
SRTp  =,  IO.56  days 
A  parametric  study  has  been  carried  out  with  r  as  parameter, 
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g  being  given  a  constant  value  of  k.O.     The  successive  values  of 

r  are  0.,  0,05,  0.10,  0.15,  0.20,  0.25  and  0.333,  respectively. 

Figure  15,  which  Is  a  plot  of  the  optimal  cost  of  treatment  versus 

the  recycle  ratio  r.   Slows  that  the  cost  is  significantly  reduced 

by  using  recycle.   Similarly,  in  Figure  16,  a  plot  of  percent  of 

reduction  of  organlcs  versus  the  recycle  ratio  r  shows  that 

recycle  increases  the  degree  of  treatment  obtained  in  an  optimum 

system.   Figure  17  which  is  a  plot  of  total  gas  production  versus 

r,  illustrates  the  increased  gas  production  obtained  by  means  of 

recycle.   Finally,  Figure  18  gives  a  plot  of  the  wash-out  liquid 

retention  tine  8  -,  as  a  function  of  the  recycle  ratio, 
"3 

As  compared  to  the  conventional  anaerobic  process,  the 
results  obtained  in  this  work  show  that  the  anaerobic  contact 
process  has  several  advantages  under  the  optimal  condition.   For 
a  recycle  ratio  of  0.25  and  a  clarifier  efficiency  of  4.0,  the 
cost  reduction  Is  (11.333  -  5.10^11.333  -  5^.9%.      In  the  conven- 
tional process  the  concentration  of  organlcs  in  the  effluent 
stream  is  0.289  +  0.386  +  1.885  =  2.56  gm/1.   The  reduction  in 
organlcs  concentration  is  then  (10.  -  2.56)/10  »  r/h,h%.      In  the 
contact  process  the  concentration  of  organlcs  in  the  outlet 
stream  is  0.289  +  0.23*1  +  0.91*  =  1.^37  ga/1  which  yields  a. 
reduction  in  organlcs  of  (10.  -  1.^3?)/10.  =  85.6;?.   The  amount 
of  gaseous  products  is,  for  the  contact  process,  (1.  +  .25)  x 
2.807  -  3.51  gm/1.   Thus,  there  is  an  improvement  in  production 
of  gaseous  px-oducts  which  can  be  used  to  meet  the  energy 
requirements  of  the  process.   This  improvement  is  (3.51  -  2.667)/ 
2.667  =  31.6%.      Finally,  the  wash-out  tiirse  ?    decreases  from 
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7.089  days  for  a  conventional  system  to  1.772  days  for  the  contact 
process. 

Figures  15  to  18  show  that  these  improvements  in  the  perform- 
ance  of  the  system  take  place  over  the  whole  permissible  range  of 
the  recycle  ratio.   However,  the  assumption  of  a  constant  value 
for  the  clarlfier  efficiency  g  may  not  be  true  when  r  takes'  on 
high  values.   At  these  recycle  rates,  almost  all  the  organisms 
contained  in  the  stream  leaving  the  second  stage  have  to  be 
separated.   Indeed,  the  major  problem  arising  from  the  use  of  the 
anaerobic  contact  process  to  date  is  related  to  an  inability  to 
separate  efficiently  the  bacterial  solids  from  the  effluent 
stream  for  recycle  back  to  the  second  stage  (14).   High  efficiency 
is  necessary  to  maintain  the  required  long  sludge  retention  time 
while  operating  at  short  hydraulic  detention  times  (14).  lis  the 
successful  full-scale  treatment  of  meat-packing  wastes  (15),  a 
vacuum  degasifier  has  been  used  between  the  digester  and  final 
settling  tank  to  remove  gases  which  tend  to  float  the  solids 
rather  than  allowing  then  to  settle  in  the  settling  tank.   This 
scheme  or  some  other  solids  separation  device  may  be  needed  to 
implement  the  anaerobic  contact  process. 

Moreover,  high  values  for  the  recycle  ratio  imply  long 
solid  retention  times.   This  yields  organism  decay,  which  is 
assumed  to  be  negligible  in  this  work.   Hence,  the  curves  given 

*  A  flotation  process  making  use  of  the  large  quantities  cf 

dissolved  gases  to  float  and  concentrate  the  solids  for  return 
to  the  digester  also  appears  feasible  (10). 
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on  Figures  15  to  18  are  not  significant  for  high  values  of  r. 
However,  they  show  a  realistic  trend  when  the  value  of  r  is 
consistent  with  the  assumptions  of  constant  clarifler  efficiency 
and  constant  maximum  specific  growth  rate. 

VI .   CONCLUSIONS 

A  mathematical  formulation  of  the  kinetic  model  reported  by 
several  authors  (7,  9,  10)  for  the  anaerobic  digestion  process 
allowed  us  to  analyze  the  two-stage  systems  under  consideration 
in  this  work.   It  has  then  been  shown  that  the  wash-out  times 
are  important  design  parameters.   In  addition,  the  process  has 
been  optimized  by  considering  an  objective  function  of  economic 
type.   Under  the  assumptions  which  have  been  made,  the  optimiza- 
tion problem  has  one  and  only  one  minimum  solution. 

As  reported  by  Pfeffer  (13) ,  this  work  shows  that  the 
contact  process  has  several  advantages  over  the  conventional 
process.   They  are  more  efficient  organic  content  reduction, 
increased  total  gas  production  and  greater  stability. 

This  work  is  a  theoretical  one.   Hence,  additional  experi- 
mental work  is  required  to  verify  the  assumptions  which  have 
been  made.   However,  since  industrial  biological  waste  treatments 
involve  mixed  cultures  metabolizing  mixed  substrates,  this  study 
constitutes  a  contribution  to  improve  our  knowledge  and  control 
of  these  processes. 
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KOMESCLATUHE 


q  Flow  ratd 


S  Concentration  of  organics  in  the  feed  stream. 

Sj         Concentration  of  primary  substrate  in  first  (and  in 
second)  stage. 

R^i  U1      Concentration  of  volatile  acids  R  and  U  in  first  stage. 

H  ,  U       Concentration  of  volatile  acids  P.  and  U  in  the 
second  stage. 

P  Concentration  of  final  product  in  the  second  stage. 

A  Concentration  of  organism  A  in  the  first  stage. 

B,  C        Concentration  of  organisms  3  and  C  in  the  second  stage. 

k , ,  kB>  k„  Maximum  specific  growth  rate  for  organisms  A,  B,  C, 
respectively. 

K  ,  Kn,  KTI  Saturation  constants  for  S,  H,  U,  respectively . 
S    n    u 


n 

3 
n  opt 

9wl 

ew2 

@w3 


Residence  time  at  stage  n. 
Optimal  residence  time  of  stage  n. 
Wash-out  retention  time  for  organ!  s:,".  A. 
Wash-out  retention  tine  for  organism  E. 
Wash-out  retention  time  for  organism  C. 


Y.,,,        Yield  constant  for  formation  of  organism  A  in  terms 


A/S 


of  utilization  of  substrate  S, 


Yo/„  Yield  constant  for  formation  of  R  from  S. 

Y..  ,„  Yield  constant  for  formation  of  U  from  S. 
J/S 

Y  .  Yield  constant  for  formation  of  organism  B  from. 


B/R 


utilization  of  intermediate  R. 


Y    ,,,  Yield  constant  for  formation  of  organism  C  from 

'  v  utilization  of  intermediate  U. 


lP/U 


Yield  constant  for  formation  of  product  P  from 
intermediate  R. 

Yield  constant  for  formation  of  product  P  from 
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intermediate  U. 

Recycle  ratio  for  the  second  stage. 

Clarify er  efficiency. 

Transformed  constants  for  kD  and  k,, . 

Transformed  yield  factor  for  YR/R- 

Transformed  yield  factor  for  Y  .  . 

Transformed  yield  factor  for  Y.,.., 

Transformed  yield  factor  for  Yp/n> 

Contact  factor. 

Rate  of  production  of  organism  A  in  the 

1  first  sts 

HjS, 

Hates  of  prcduction  of  organisms  B  and 
in  the  second  stage. 

C  respect J 

voly 

Rate  of  consumption  of  substrate  S  in  1 

:he  first-  ,< 

jtage, 

Rates  of  consumption  of  organics  R  and 
in  the  first  stage. 

U  respect! 

vely 

Rates  of  consumption  of  organics  R  and 
in  the  second  stage. 

U  respeptl 

.vely 

Rate  of  production  of  product  P  in  the 

seccnd  ct:: 

IgjB. 

Objective  function. 

Penealization  for  discharging  organics 
respectively. 

S ,  H ,  and 

u 

Response  surface. 

Minimum  point  on  (I). 

Saddle  points  on  (E). 

Maximum  point  on  (li). 

Kessian  matrix. 
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COMPUTER  PROGRAM 


Main  routine  MICU20   and   subroutines 
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C  MAIN  PROGRAM  MICU20. 

C  THIS  PROGRAM  DEFINES  THE  OPTIMAL  DESIGN  OF  A  TWO-STAGE  WATER  TREATMENT 

C  BY  MIXED  CULTURES.  IT  CALLS  k   SUBROUTINES  SIMPLE, OBJ .NEWTON, IS 03AT 

C  ■  1  FUNCTION     QUART 

C  PROGRAMMER   JC  DALTES  KSU, MANHATTAN, FALL  63 

C 

C  FIRST  PART*«***AFECIFICATI0N3,  STATEMENT  FUNCTION  DEFINITIONS  AND  DATA 

C  READING . 

EXTERNAL  03J 

REAL  KA,KS,KB,KR,KC,KU,MUS,MUS,MUU,KBO,KCO 

DIMENSION  TITLE(30) ,TH(2,3) ,TMIN(2) ,DV(2) 

COMMON   TH171,THW2,THW3,KUS  ,MUR,MUU,KA,KS  ,KB,KC  .KU.DELTS  .DELTAY 

COKMON   /SEARCH/MDIK ,KDIM ,TH .ALPHA, BETA , GAMMA , IBARY , EP3IL , LIMIT , 
1  ICONV.TMIN.SMIN, NITER, CHIT.NEVAL 

C 
C    READ   AMD   ECHO-CHECK   DATA.    DEFINE   STATEMENT    FUNCTIONS. 

XlSl(T)    =   KS/(T*KA   -   1.) 

DELTAS (T)=  X1SO  -  X1S1(T) 

X2A1 ( T )=  YAS*DELTAS ( T ) 

X1RKt)  =   YRS*DELTAS(T) 

X1U1(T)=   YUS*DELTAS(T) 

X1R2(T)=   KR/(T»KB   -    1.) 

DELTAR(T,U)=   X1R1(T)-X1R2(U) 

X1U2(U)   =   KU/(U*KC-1.) 

DELTAU(T.U)   =  XIUI(T)    -  X1U2(U) 

X2B2(T,U)   =  YBR*DELTAR(T.U) 

X2C2(T,U)   =  YCU*DSLTAU(T,U) 

X1P2(T,U)   =   YPR/YBR*X?.B2(T,U)   -:■   YPU/YCU*X2C2(T,U) 

W2(T1)=    (KR+X1R1(TM/(KB*X1R1(T1) ) 

W3(T1)=,    (KU+X1U1CT1)  )/(KC*XlUl(Tl) ) 

REAE( 1,1010)      TITLE 

WRITE (3, 3010)    TITLE 

READ(1,1020)    X1S0,KA,KS,YAS,YRS,YUS 

WRITE(3,3020)X1SO,KA,KS,YAS,YRS,YUS 

READ( 1,1020)    KS0,KR,Y3R0,YPR0 

WRITE(3,3030)K3O,KR,YBRO,YPRO 

EEAD( 1,1020)  KCO.KU.YCUO.YPUO 

WRITE ( 3 , 30^0 ) KC  0 ,KU , YCUO , YPUO 

RE AD ( 1 , 1 0  2  0 )      MUS , MUR , MUU 

WRITE(3i30ol)  MUS, MUR, MUU 

PARAK=0 
5   READ(1,1020,END=S5)    R.BET 

CONTAC=   l.-lBbT-i. J*R 

WRITE(3,3075)    R.BET.CONTAC 

IF(CORTAC.GT.O.O)    GO  TO  10 

WRITE(3,30',;6) 

STOP 

r 

C   GENERATE   FICTITIOUS   CONSTANTS   FOR  THE  SECOND  STAGE. 
10  KBbKBO/CONTAC 
K&dCCO/COHTAC 
YBR.-YBHO/CCNTAC 
STCOWfCUO/CONTAC 


66 


YPIfcYPRO/a.-j-R) 

YPU=YPUO/(l.+R) 

WHITE(3,3077)  K3,KC,YBR,YCU,YPR,YPU 

THW1=  ( KS+X13  0 )  /  ( KA*X1S  0 } 

IF(FARAM.EQ.l)  GO  TO  12 

IF(NOFT.EQ.O)  GO  TO  25 
C 
C  IF  OPTIMIZATION  PROBLEM,  READ  IN  AND  WRITE  PARAMETERS  FOR  SEARCH. 

WRITE(3,3083) 

NDIK=2 

KDIMbNKOW. 

WRITE(3i3050)NDIM 

READ( 1 , 1070 )  ALPHA , DETA, GAMMA , STEP , IBARX 

WRITE (3, 3061) 

WRITE (3,3070) ALPHA , BETA , GAMMA , STEP , IBARY 

READ(1,1080)    EPSIL, LIMIT, ETA, MAXIT 

WRITE(3,30S0)    E3PSIL.ETA 

WRITE (3, 308 2) 
C 

C   SECOND   PART*****DELTERMINATION   OF   OPTIMAL  POLICY 
C   DIFFERENTIAL   CALCULUS    APPROACH. 
12  WRITE(3,3085) 

TH10PIW1 . /KA* ( 1 .+SQRT ( MUS*KS*KA) ) 

TBSTls   MUS<<KA*X1S0*X1SG  -   KS 

CAl,L  NEWTON (TH20PT) 

THV73    (KUvXlUl(THlOPT)  )/(KC*XlUl(TH10PT) ) 

TEST2=    QUART (THW3) 

WRITE ( 3. 3086 )    TH10PT.TEST1 .TE20PT ,TEST2 

IF(TEST1.GT.0.0.AND.TE3T2.GT.0.0)    GO  TO   15 

WRITE (3, 3088) 

STOP 
C 

C   SEARCH  TECHNIQUE  APPROACH. 
15  WRITE (3, 303 7) 

LLAC=0 

NEVAL=0 

HTH1=  TH10PT+0.5 

BTHZa  TH20PT+2 . 
C 

C    ENTER   OPTIMIZATION   LOOP. 
21   TH(1,1)=   TfilOPT+STEP 

TK(2,1)=  TH20PT 

TH(1,2)=  TH10PT 

TH.'2,2)  =  TH20FT+STEP 

TH ( 1 , 3 )=  TH10PT+STEP 

TH(2,3)=  TH20PT+STEP 

THW2=  W2(BTH1) 

THW>_-  W3(P.TH1) 

CALL  SIKPLE(OBJ) 

LLAC=I,LAC+1 

IF(ICONV.EQ.l)    GO  TO   22 

WRITU'(3,3150)    NITER, CHIT 

STOP 
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22  IF(A3S(TMIN(1)-RTH1).LT.ETA.AND.ABS(TMIN(2)-RTH2).LT.ETA)  GO  TO  23 
IF(LLAC.EQ.MAXIT)  GO  TO  2k 

RTH1-=TKIN(1) 
RTH2=THIN(2) 
GO  TO  21 
2*  WRITE ( 1,3100)  LLAC 
STOP 

23  BTla  TMIN(l) 
HT2=  TMIN(2) 
THW2=  W2(RT1) 
THW3=  W3(RT1) 
WRITE ( 3.3109)  LLAC 
WRITE ( 3 , 3110 )  RT1 , RT2 
WRITE (3, 31^0)  SMIN 
GO  TO  26 

C 

C  THIRD  PABT*****THE  POLICY  IS  FIXED  BY  THE  USER  WITHOUT  BEING  OPTIMAL. 
C  IP  MO  OPTIMIZATION  PROBLEM 
2  5  WRITE(3,308i|-) 

READ(1,1020,END=999)  RT1.RT2 

WRITE (3, 30*1-2)  BT1.RT2 

DV(1)=  RT1 

DV(2)=  RT2 

THW2=  W2(RT1) 

THW3=  W3(RT1) 

NEVAI..0 

WRITE (3. 308 2) 

CALL  0&J(DV,S) 

WRITE(3.3;4-90)  S 
C 

C  FOURTH  PART*****  STATE  CORRESPONDING  TO  THE  CHOSEN  POLICY 
26  WRITE(3,3500)  THW1.RT1 

IF(RTl.GT.THWl)  GO  TO  Jk 

30  WRITE(3,3079) 

STOP 
3*+  Sl^  XlSl(RTl) 

Al=  X2A1(RTI) 

Rl=  XlRl(RTl) 

Ul=  XlUX(RTl) 

WRITE(3,3120)  Sl.Al.Rl.Ul 

WRITE ( 3 , 3520 )  THW2 ,THE3 , RT2 

IF(RT2.GT.THW3)  GO  TO  35 

GO  TO  30 
35  R2=  X1R2(RT2) 

U2=  X1U2(RT2) 

B2=  X2B2(RT1,RT2) 

C2=  X2C2(RT1,RT2) 

P2=  X1F2(RT1,RT2) 

WRITE(3,3130)  R2,U2,B2,C2,P2 

SRT-_=    RT2/CONTAC 

WBITE(3,3089)    SHT 

IF(ICUBV.EQ.O)    GO  TO   75 
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C   FIFTH  PART  •**•*   OUTFUT   OF  CURVES   CORRESPONDING  TO  THE  CHOSEN  FOLICY, 
C  1)    FIBST  STAGE 

36  WHITE (3, 3180) 

Slo  XlSl(THWl) 

Al=  X2A1(THW1) 

Rl=  XTRl(THWl) 

Die   XlUl(THWl) 

DX1S1=   -ICS*KA/(KA*THW1-1.  )**2 

DX2A1=    -YAS*DX1S1 

DX1R1=    -YRS*DX1S1 

DX1U.U-   -YUS*DX131 

WRITE(3,3190)    THva,Sl,Al,Rl,Ul,DX131,DX2Al,DXlRl,DXlUl 

RT  =    (AINT(10.*THW1)+1.)/10 
40   Slo   XlSl(RT) 

Al=   X2A1(RT) 

Rl=   XlRl(RT) 

Ul=   XlUl(RT) 

WRITE(3,3190)    RT,S1,A1,R1,U1 

RT=   HX+0.1 

DJJo   RT   ~  AINT(RT) 

IF(. NOT. (1.-DIF.LT.0. Ol.OR.DIF.LT. 0.01))    GO  TO  40 
50  Sla  XlSl(RT) 

Al=   X2A1(RT) 

Rl=   XlRl(RT) 

Ulx.   XlUl(RT) 

WRITE(  3 ,  3190 )    RT ,  SI ,  11 ,  Rl ,  Ul 

RT=-HT+1. 

IF(RT.LE.20.0)  GO  TO  50 
C 
C      2)  SECOND  STAGE  BETWEEN  THW2  AND  THW3. 

WRITE (3, 3200) 

R2=  X1R2(TKW2) 

B2=  X2B2(RT1,TKW2) 

P2=  YPR/YBR*B2 

DX1R2=  -KH*KB/ ( THW2*KB-1 . ) **  2 

DX2B2=  -YBR--DX1R2 

D1X1P2=  YPR/YBR*DX2B2 

WRITE( 3 , 3210 )  THW2 , R2 , B2 , P2 . DX1R2 ,DX2B2 ,D1X1P2 

RT=  (AINT(10.*THW2)+1.)/10 
60  R2=  X1R2(RT) 

B2=  X2B2(BT1,RT) 

F2=  XPB/XBB»82 

WHITE(3,32"10)  RT,R2,B2,P2 

RT=RT+0.1 

DIFrt  RT  -  AINT(RT) 

IF(. NOT. (1.-DIF.LT.0. 01.DR.DIF.LT. 0.01))  GOTO  60 
6i  H2^  X1R2(RT) 

B2=  X2B2(RT1,RT) 

P2=  YrR/YBR*B2 

WRITE(3, 3210)  RT,R2,B2,P2 

BI=HD+1. 

IF(RT.GE.TRW3)  GO  TO  ?1 

GC  TC  61 
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C  3)    SECOND  STAGE  ABOVE  THW3 

71  WHITE(3,3220) 
R2=  X1R2(THW3) 
B2=  X2B2(RT1,THW3) 
P2=  X1P2(RT1,THW3) 
U2=  X1U2(THW3) 
C2=  X2C2(RT1,THW3) 

EKlUa*   -KU*KC/(THVJ3*KC-1.  )**2 
DK2C2=   -XCU*EX1U2 

D2XlP2-=   YPR*KR*KC/(KB*THW3-1.  )**2 
D3X1P2=   YPU/YCU*DX2C  2+D2X1P2 

WRITE(3,3230)    THV.'3IR2,B21P2,U2,C2,DX1U2,DX2C2,D2X1P2,D3X1P2 
BT=    (AIKT(10.*THW3)+1.)/10 
73  B2=  X1R2(RT) 

B2r=   X2B2(RT1,RT) 

P2=  X1P2(RT1,RT) 

U2=  X1U2(RT) 

C2=   X2C2(RT1,RT) 

WRITE(3,3230)    RT,R2,B2,P2,U2,C2 

KD«rHT+G.l 

DIF=:RT-AINT(RT) 

IF( .NOT. (l.-DIF.LT. O.Ol. OR. DIF.LT. 0.01))    GO  TO  73 

72  R2=  X1R2(RT) 
B2=  X2B2(RT1,RT) 
P2_   X1P2(R'TJ.,RT) 
U2=  X1U2(RT) 
C2=  X2C2(RT1,RT) 

WRITE(3,3230)    RT,R2,B2,P2,U2,C2 
BT=BT+1. 

IF(RT.LE.20.)    GO  TO   72 
75   IF(ICONT.NE.O)    GO   TO   80 
PABAM=1 
GO  TO   5 
C 

C    SIXTH   PART   »«*»•    DRAW  THE  CONTOURS    OF  THE   OBJECTIVE   FUNCTION. 
80   HEAD(1,1020)    DELTS.DELTAY 

CALL   IS03AT 
85  WRITF.(  3,32^0) 
GO  TO  1000 
999   WR1TE(3,35^0) 
1C00   STOP 

1010    FORMAT (  20A^/10A'0 
1020    FOBMAT(6F10.0) 
1041   FORMAT (3110) 
1070   FORMAT  ( if  -PIC  .0,110) 
1080    FORMAT(2(E10.2,110)) 

3010   FORMAT ( 110 , 30AV///T55 , 5H*****/T55 , 6H*DATA*/T55 t 5H***»*//) 
3020    FORHAT(19H   FEED  CONCENTRATION  ,  21X,'+KX1S0,F9. 2, 5H   GK/L/ 

1  29H   MAXIMUM  SPECIFIC   GROWTH   RATE,11X,2HKA, F11.2 ,6H  DAY-1/ 

2  20H   SATURATION   CONSTANT, 20X.2HKS .F11.2.5H   GM/L/ 

3  l^H   YIELD   FACTORS, 26X,3HYAS,F10. 2/ 
*t  40X,3HYR3,F10.2/ 


?0 


5  40X.3HYUS.F10. 2//) 

30.30    FORMAT (29H   MAXIMUM  SPECIFIC    GROWTH   RATE ,  11X ,  2HKB ,  Fll .  2 ,  6H   DAY-1/ 

1  20H   SATURATION   CONSTANT , 20X.2HXR, Fll. 2, 5H   GM/L/ 

2  14H  YIELD   FACTORS,  26X,3HYBR,F10. 2/ 

3  40X,3HYPR,F10.2//) 

3040    FORMAT (29H   MAXIMUM   SPECIFIC    GROWTH   RATE ,  11X  ,  2HKC  ,  Fll  .  2 ,  6H   DAY-1/ 

1  20H  SATURATION  CONSTANT.  20X.2HKU, Fll. 2,  5H  GM/L/ 

2  14H  YIELD   FACTORS, 26X,3HYCU,F10. 2/ 

3  40X.3HYPU.F10. 2//) 

3042  FORMAT (16H  RESIDENCE  TIMES , 24X.3HRT1 .F10.2/ 

1       40X.3HRT2.F10. 2/) 
3050  FORMAT (18H  SEARCH  PARAMETERS/15H0DIMEN3IONALITY ,25X,4HMDIM,l6/) 
3061  FORMAT (9H  STRATEGY/) 
30?0  FORMAT (11H  REFLECTION, 29X, 5HALPHA.F8 .2/ 

1  12H  CONTRACTION  28X.4HBETA, F9. 2/ 

2  10H  EXPANXION  30X.5HGAMMA.F8 . 2/ 

3  T4l,'STEP' .F9.2/T41.5HIBARY.15/) 

3075  FORMAT ( '1RECYCLE  RATIO' ,T4l, 'R' , F13.3/ 

1  'CLARIFYEH  EFFICIENCY' ,T4l, 'BET' , Fll. 3/ 

2  'CONTACT   FACTOR'  ,T4l ,  "CONTAC  '  ,  F8  . 3//) 

3076  F0RMAT('1THE   CONTACT   FACTOR   IS   NOT   POSITIVE.'//) 

3077  FORKAT( 'FICTITIOUS   CONSTANTS  1    KB/KC/YBR/YCU/YFR/YPU' ,6f10. 2//) 

3079  FORMAT( ' 0THI3   POLICY    IS    NOT   FEASIBLE'///) 

3080  F0RMAT(20H   PRESCRIBED   ACCURACY ,20X ,5HEPSIL,lPEll.l/ 

1  33H   ACCURACY    OF   OVER   ALL  CONVERGENCE, ?X, 3HETA.1PE13.1/) 

3061   FOHMAIOIH  PENALIZATION   FOR  DISCHARGING  S , 9X , 3HMOS , P12 . k/ 

1  30X,1HR,9X,3HHUR,F12.4/ 

2  30X .1HU.9X , 3HMUU.F12.4//) 

3082  FORMAT ( 1K1  ,T56  ,  9H****«****/T56 ,  9HttRESULTS*/T56  ,  9H***»*«***/// 
1  47H  CONCENTRATIONS   ARE   IN  GM/L.    TIMES    ARE   IN  DAYS.///) 

3083  FORMAT(36HODETSRMINATIOH   OF  THE   OPTIMAL  POLICY///) 

3084  FORMAT (33HOTHE  CHOSEN   POLICY    IS    NOT    OPTIMAL/ 

3085  FORMAT (46X.6H*****   14HOPTIMAL  POLICY   6H  *•***// 
1  45X.30HDIFFERENTIAL  CALCULUS    APPROACH//) 

3086  P0RMAT(1H    39X.6HTK10PT.1PE12.3/1H    39X, 5HTEST1 .1PE12. 2// 
1  1H   39X,6KTH20PT,1PE12.3/40X,5HTEST2,1PE12.2//) 

308?    FORI i AT (5 OX,  16HS IMPLEX   APPROACH//) 

3033    FORMAT ( 'INEGATIVE   TEST    FOR   FEASIBILITY.    THE   CALCULATIONS    ARE', 

1  *    STOPPED'//) 

3039    FORMAT ( /// '    SOLID   RETENTION   TIME    IN   SECOND   STAGE', 

1  T4l,'SRT',F11.3//) 

3100   FORMAT(31H  CALCULATIONS   ARE  STOPPED  AFTER, 15, 12H   ITERATIONS.) 
3109   POHMAT(35H  CONVERGENCE  HAS   BEEN   REACHED  AFTER, 14, 11K   ITERATIONS/) 
3116    FC2MAT(27H    OPTIMAL   DFCISICN   VARIABLES  .1PE3I.3/1H   1PE57.3/) 

1  40X ,  4HX2A1 ,  F10  .  3/4 OX ,  4,1X1111 ,  PIC .  3/40X ,  4FX1U1 ,  F.I  0 .  3/ ) 

3130    FORMAT  (33:1   STATE   VECTOR   AT   THE   SECOND   STAGE, 7X,4KXlR?.Flb.  3/ 

1  40X , 4HX1U2 , F10 . 3/40X , 4HX2B2 , F10 . 3/40X , 4HX2C  2 , F10 . 3/ 

2  40X.4HX1P2.F10.3/) 

3140   FORKAT(2?H   OPTIMAL   OBJECTIVE   FUNCTION, 1PE32. 4/) 

3150    FORKAT(6ti    AFTER, 15, 33H    ITERATIONS,    THE  CRITERION   EQUALS    1PE10.2/ 
1      2bh  CALCULATIONS   ARE  STOPPED./) 
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3180  FORMAT ( Utt, ^3X,6H****«  20HCONCENTRATION  CURVES  6H  »♦**•///// 

1        55X.11HFIRST  STAGE/// 

1  5X ,  2HTH ,  10X ,  1J-HX1S1 ,  11X ,  4HX2A1 ,  11X ,  *}HX1R1 , 

2  3 IX , teXlUl , 10X , 5HDX1S1 , 10X , 5HDX2A1 , 10X , 5HDX1H1 ,10X, 5HDX1U1/ } 
3190  FOHKAT(F8.3,F13.3i7F15.3) 

3200  FORMAT (lHl,'!-5X,3'ffiSECOND  STAGE  BETWEEN  THH2  AND  THW3//// 

1  kX ,  2HTH  ,  10X ,  'I-HX1R2 , 

2  11X , kSXZBZ ,11X, ^HX1P2 , 10X , 5HDX1H2 , 10X , 5HDX2B2 , 10X , 6HD1X1P2/ ) 
3210    F0RKAT(F?.3.F1'I-.3.5F15.3) 

3220   FOBMAT(1H1,50X,23H3ECOND  STAGE  ABOVE  THW3//// 

1  bX. ,  2RTH ,  ]  OX  ,  4HX1R2 ,  11X ,  taX2B2 , 

2  llX,4HXlP2,llX,'+HXlU2,12XltaX2C2,    7X.5HDX1U2,    7X.SHDX2C2,    6X, 

3  6HD2X]F2,6X,6HD3X1P2//) 

3230  FORMAT (F7.  3.  F14.  3,^15.3,^12. 3) 

32^0  FORMAT (19H1N0RHAL  TERMINATION) 

3^90  FORMAT (19H  OBJECTIVE  FUNCTION, F35. 3/) 

3500  F0RMAT(1H0,52X,11HFIRST  STAGE//14H  WASH  OUT  TIME.F'MJ .3// 

1   15H  RESIDENCE  TIME.F39.3/) 
3520  FORMAT(lHO,52X12HSECOND  STAGE//15H  WASH  OUT  TIMES.F39.3/ 

1        1H  F53.3//15H  RESIDENCE  TIME.F39.3/) 
35^0  FORMAT (  31H1NO  RESIDENCE  TIMES  FOUND.  STOP) 
END 
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SU3R0UTINE  S IKPLE ( OBJ ) 
C  THIS  SUBROUTINE  FINDS  THE  MINIMUM  OF  THE  OBJECTIVE  FUNCTION  GIVEN 
C  BY  SUBROUTINE  OBJ.  THE  ARGUMENT  OF  SIMPLE  REQUIRES  AN  EXTERNAL 
C  STATEMENT  IN  THE  CALLING  ROUTINE. 

C  FOLLOWING. VARIABLES  SHOULD  BE  PREVIOUSLY  DEFINED  AND  STORED  IN 
C  COMMON  /SEARCH/ 

C  NDIM.KDIM , TH , ALPHA , BETA , GAMMA , IBARY , EPSIL , LIMIT 
C  THIS  S.NE  DEFINES  AND  STORES  INTO  THE  SAME  AREA 

C  ICONV.TMIN.SMIN, NITER, CHIT "WHERE  THIN  IS  THE  OPTIMAL  DECISION  BECTOR. 
C  NEVAL  IS  DEFINED  BY  S.NE  OBJ. 

C  IP  THE  SEARCH  CONVERGED,  ICONV=l.  IF  NOT,  ICONV=-l  AND  THE  BEST 
C  CURRENT  POINT  IS  CONSIDERED  AS  OPTIMAL. 

C  DIMENSIONS    T,SIGTH,TBAR,TREF,TEXP,TCON ,TMIN      NDIM 
C  S,W  KDIM 

C  TH  NDIM.KDIM 

C  CHANGING  THE  DIMENSIONALITY  OF  THE  SEARCH  REQUIRES  THE  PROPER 
C  DIMENSION  STATEMENT.  THIS  SUBROUTINE  HAS  BEEN  TESTED  FOR  NDIM=1 , 2, 3,4, 20. 
C  PROGRAMMER  JC  BALTES  KSU,  CHEM  ENG,  FEB  68' 

C 
C  SPECIFICATIONS. 

COMMON  /SEARCH/NDIM ,KDIM , TH , ALPHA , BETA , GAMMA , IBARY , EPS IL , LIMIT , 
1  IC  ONV , TMIN , SMIN , NITER , C  HIT , NEVAL 

DIMENSION  TH(2,3) ,T(2) ,S(3) ,SIGTH(2) ,TBAR(2) ,TREE(2) ,TEXF(2) , 
1  TCON(2) ,TMIN(2),W(3) 

C 
C  INITIALIZATION 

NITEB=0 

ICONV=0 
C 
C  COMPUTE  INITIAL  FUNCTION  VALUES. 

DO  30  J-l.KDIM 

DO  40  1=1, NDIM 

40  T(1)=TH(I,J) 

30  CALL  03J(T,S(J)) 
C 

C  BEGINNING  OF  ITERATIONS. 
C  DEFINE  THE  POINT  HAVING  THE  HIGHEST  FUNCTION  VALUE 

41  SH  -.:     S(l) 

JHhI 

DO  50  J=2,KDIM 
IF(S(J).LE.SH)  GO  TO  50 
SH=S{J) 
JHbJ 

50  CONTINUE 
C 

C  DEFINE  THE  POINT  HAVING  THE  LOWEST  VALUE  OF  THE  OBJECTIVE  FUNCTION. 

51  SL=S(1) 

JIc:l 

DO  60  J^2,KDIM 
IF(S(J).GS.SL)  GO  TO  60 
SI«S£J) 

JL-rJ 

60  CONTINUE 
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IF(IC0NV*IC0NV.EQ,1)  GO  TO  320 
C 

C  DETERMINATION  OF  THE  CEKTROID 
IF(IBARY.NE.O)  GO  TO  90 
DO  70  Ir.-a.NDIM 
SIGTH(I)=0. 
DO  80  Jt=1,KDIM 
IF(J.EO.JH)  GO  TO  80 
SIGTH( I )aSIGTH(I )+TH( I , J ) 
80  CONTINUE 
70  TBAR(l)=rSIGTH(I)/FLOAT(NDIM) 

GO  TO  130 
90  SIGMAlfeO. 

DO  100  J=1,KDIM 
W(J)=S(JH)-S(J) 
100  SIGMAW=SIGMAW+W( J ) 
DO  110  I=1,NDIM 
SIGTH(I)=0. 
DO  120  J=1,KDIM 
120  SIGTH(  I  )=SIGTH(  I  )+W(  J  )*TE(  I. ,  J  ) 
110  TEAR ( I  )--SIGTH(  I  )/SIGMAW 
130  CALL  OEJ ( TEAR, S BAR) 
C 

C  REFLECTION  GIVES  POINT  (TREF) 
DO  140  1=1,NDIM 
I'lO  TBEF(  I)=  (1 .+ ALPHA) *T3AH(  I )  -AlPHA*TH(I,JH) 
CALL  OBJ (TBEP, SEEP) 
IF  (SRHF.LT.S(JL))  GO  TO  160  . 
DO  150  J=1,KDIM 
IF(J.EQ.JH)  GO  TO  150 
IF  (STREF.GT.S(J))  GO  TO  150 
GO  TO  180 
150  CONTINUE 
GO  TO  220 
C 

C  EXPANSION  GIVES  POINT  (TEXP) 
160  DO  170  I=1,NDIM 

170  TEXP(I)=  GAKI!A*TREF(I)  +  (l.-GAMMA)*TBAR(I) 
CALL  OBJ ( TEXP , SEXP ) 
IF(SEXP.LT.S(JL))  GO  TO  200 
180  DO  190  I=1,NDIK 
190  TH(I,JH)=TREF(I) 
S(JH)=S3EF 
GO  TO  300 
200  DO  210  I=1,NDIM 
210  TH(I,JH)=TEX?(I) 
S(JH)=SEXP 
GO  TO300 
C 

C  CONTRACTION  GIVES  POINT  (TOOK) 
220  IP<8HEP.GT.S(JH))  GO  TO  2'+0 

DO  230  Ial.NDIM 
230  toi{I,JH}=TBEP(I) 


7* 


2^0  DO  250  Ir.-1,NDIM 

250  TC0N(1')=3ETA*TH(I,JH)  +  (1.-  BETA)*T3AR(I) 

CALL  OBJ (TC ON, SCON) 

IF(SCON.GT.S(JH) )  GO  TO  270 

DO  260  I=1,NDIM 
260  TH(I,JH)=TCON(I) 

S(JH)=SCON 

GO  TO  300 
2?0  DO  290  J=1,KDIM 

IF(J.EQ.JL)  GO  TO  290 

DO  250  I=1,NDIM 

TB(I,J)=  (TH(I,J)+TH(I,JL)0/FLOAT(NDIM) 
280  T(I)=TH(I,J) 

CALL  OBJ(T,S(J)) 
290  CONTINUE 
C 

C  CRITERION  FOR  CONVERGENCE. 
300  SQDELS=0. 

DO  310  J=1,KDIM 
310   SQDELS=SQDELS+(S(J)-S3AH)**2 

CRIT=SQRT(SQDELS/FLOAT(NDIM) 

NITER=NITER+1 

IF(CRIT.GS.EPSIL)    GO   103*4-0 

ICGSVatl 

GO  TO   51 
3^0    IFiKITER.LT. LIMIT)    GO   TO   kl 

ICONV=-l 

WRITE (3,3010)  NITER 
3010  FORMAT C+9H1 THE  SEARCH  RETURNS  THE  BEST  CURRENT  POINT  AFTER, 

1        12HITERATION  NO  15.31H.  THIS  IS  NOT  THE  TRUE  OPTIMUM.//) 

GO  TO  51 
320  DO  330  I=1,NDIM 
330  THIN(I)=TH(I,JL) 

SMIN=S(JL) 

RETURN 

END 
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FUNCTION  QUAR(U) 

REAL  KA,KS,K3,KR,KC , KU , MUS , MUR , MUU 

COW.HON  THW1  ,THW2,THW3,MUS  ,MUR,MUU,KA,KS  ,KB,KR,KC  .KU.DELTS  .DELTAY 

QUART       =     MUB*KB«KB/{U*KB-1.}*»2  +   MUU*KU*KC/(U*KC-.1«  0**2-1. 

RETURN 
END 


76 

SUBROUTINE  NEWTON ( ROOT ) 
C 
C  THIS  SUBROUTINE  COMPUTES  A  ROOT  OP  QUART=0 

REAL  KA , KS , KB , KR , KC , KU , KUS , HUR , MUU 

COMMON  THW1,THW2,THW3,MUS ,KUR,MUU,KA,KS .KB.KR.KC .KU.DELTS .DELTAY 

DQUART(U)=  -2.*MUR*KH»KB*KB/(U*KB-1.)**3 
1        -2.*MUU*KU*KC**2/(U*KC-1.0**3 
C 
C  LOOK  FOR  AN  INTERVAL  WHERE  FUNCTION  QUART  HAS  OPPOSITE  SIGNS. 

Xl-_=  1./KC+O.01 
5  X2=  Xl+0.10 

IF(QUAIiT(Xl)*QUART(X2)  .LE.O.O)  GO  TO  10 

X1=X2 

GO  TO  5 
C 

10  X2=  XI  ~(QUART(X1)/DQUART(X1)) 

IP(ABS(X1-X2).LT.1.0E-0^)  GO  TO  20 

X1=X2 

GO  TO  10 
20  ROQ&X2 

RETURN 

END 
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SUBROUTINE  0BJ(T,3) 

REAL  KA , KS , KB , KH , KC , KU , MUS , MUB , MUU 

C  OKKON  THW1 ,  THW2 ,  TKW3 ,  MUS  ,  MUB ,  MUU ,  K A ,  KS  ,  KB ,  KR ,  KC  ,  KU ,  DELTS  ,  DELTAY 

COMMON   /SEARCH/   BICON( 20) .NEVAL 

DIMENSION   T(l) 

IF(T(1).LT.THW1.0R.T(2) .LT.THW3)    GO  TO  10 

X1S2=  KS/IKA*T(1)-1.) 

X1B2=  KR/(KB*T(2)-1.) 

X1U2=   KU/ ( KC«T ( 2 ) -1 . ) 

Sm  (T(l)+T(2))   +   MtJS*X182  +   MUR*X1R2  +   MUU*X1U2 

NEVAIc   NEVAL+1 

GO  TO  20 
10   S-.   l.OE+06 
20   RETURN 

END 


SUBROUTINE  ISOBAT 
C 
C  SPECIFICATIONS. 

HEAL  KA,KS,KB,KE.KC,KU,HUS,MUR,MUU 

INTEGER  SilO 

DIMENSION  TITLE(30) ,TH(2,3) ,TMIN(2) ,DV(2) 

COMMON   IHW1 ,  THVJ2 ,THW3 , MUS  , MUR , MUU ,KA ,KS  ,KB ,KH ,KC  ,KU .DELTS  .DELTAY 

C OHMON   /SEARCH/HDIM , KDIi-I , TH ,  ALPHA , BETA , GAMMA , IBARY  , EPS IL .LIMIT , 
1  ICONV ,TMIN , SHIN , NITER ,CRIT , NEVAL 

C 

WHITE(3,32ilO) 

S=      (AIKT(30.*SMIN)+1.)/10. 

S110=0 
76  DO  110   1=1,10 

WRITE(3,3250)    5 
C 
C    INCREASE  Y. 

Y=    (AINT(TIIIN(2)*10.)+1.)/10 

DO  80  J=l,25 

IF(A3S(Y-l./KB).LT.l.OE-04.OR.ABS(Y-l./KC).LT.l.OE-0i|-)    GO  TO   79 

A=   Y+MUR«-KR/(K3*Y-1.)    +   KUU*KU/(KC*Y-1 .  ) 

DISCR.-.-    ((A--3}*KA-1.)**2   -4.*KA*(MUS*KS.!.S-A) 

1F(DISCR.LT.0.)    GO  TO  81 

XX-,    (l.-(A-S)*KA  +   SQRT(DISCR))/(2.»KA) 

X2=    (l.-(A-3)*KA   -   SQRT(DISCR))/(2.«KA) 

WHITE( 3, 3260 )X1  ,X2,Y 

79  XWY+DELTAX 

80  CONTINUE  ' 
C 

C    DECREASE  Y. 

81  Y=    (AINT(TKIN(2)*10.))/10. 
DO  90   J=l,20 

IF(ABS(Y-l./KB).LT.l.OE-0iKOR.ABS*Y-l./KC).LT.l.OE-0(f-)    GO  TO  89 
A=  Y+KUR*KR/(KB*Y-1.)   +   MUU*KU/(KC*Y-1.) 

DISCR=    ((A-S)«KA-1.)**2   -^.^KA«-(HUS^KS+S-A) 
IF(DISCR.LT.O.)    GO  TO  100 
XL=    (l.-(A-S)»KA  +    SQRT(DISCR))/(2.»KA) 
X2=    (l.-(A-S)*KA   -   SORT(DISCR)  )/(2.*KA) 
SrffiITE(3,3260)    X1.X2.Y 

89  Y=Y-DELTAY 

90  CONTINUE 
100   SaS+DSXTS 

IF(SllO.EQ.l)    RETURN 
110   CONTINUE 
S110=l 
S=1.10*SHIN 
WRITE(3,325D 
GO  TO   76 

3240  p0bmat(1h1,^9x,6h*****  8hcontouhs  6h  **•**//) 

3250  format ( 19hooeiective  function ,  fl  2 .  ?. ,  a'x  ,  2hx1 ,8x,  2hx2 , 8x ,  1hy/ ) 

3251  pcr;:at(  34hocontour  corresponding  to  i.i*shih) 
3260  format (1h  .f38.3.2f10.3) 

END 
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APPENDIX  I 


SLOPE  OF  THE  CONCENTRATION  CURVES  AT  THE 
POINTS  OF  DISCONTINUITY 


At  the  points  of  discontinuity  of  the  concentration  curves, 

wash-out  occurs  for  organism  A,  B  or  C.   In  order  to  draw  these 

curves  with  accuracy,  it  is  useful  to  determine  their  slopes  at 

the  wash-out  residence  times  6  , ,  9   ,  and  6  -. 

wj.   w2       w3 

In  the  first  stage,  wash-out  occurs  when  6=0   .   Then. 

1    wl 

from  equations  (12),  (14),  (16)  and  (17),  we  obtain  respectively 


as. 


deiKi 


"  KS  kA 


(kA  9W1  -  l> 


dA 
do! 


"wl 


fA/S  dG- 


'wl 


dR 
__1 

5et 


dS, 


3wl  =  "  R/3  d9l 


wl 


dU 


as. 


de"  e  ,     u/s  d@ " 

1 1  wl  1 


'wl 


In  the  second  stage,  wash-out  occurs  for  organism  B  and  C 


when  0  <  8   <  G 

2  -  v 


When  G   < 
w2 


!  <_  6   ,  organism  B  alone  grows 


using  organics  B  and  yielding  product  P.   Then,  from  equations 
(19)  i  (?-4),  and  (27),  we  obtain  respectively 


dEg 


K  k 

R  B 


(k   G  „ 
B  w2 


1) 
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dB 
d6„ 


'w2 


dR, 
LB/R  dj~|0 


w2 


and 


dP 


"F/H    dB 


d32|9W2"  YB/R  "  d°2l6w2 


_dp 
de„ 


P/R  dB 
°w3-=  5B/R  d92 


w3- 


When  e  >  0   ,  both  organisms  B  and  C  grow.   Concentration  U 
2    w3  2 

falls,  yielding  more  product  P.   From  equations  (20),  (25),  and 
(27),  we  obtain  respectively 


dU, 


Kukc 


w3 


(kc  ew3  "  1} 


dC_ 
d9_ 


Y 


"w3+ 


C/U  d6„ 


»3 


dP         _P/R   _dB_ 

"av*"  VR'de2 


YP/U    dC 
°W3   YC/U  "  dS2  [ew3 
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APPENDIX  II 


GRAPHICAL  DETERMINATION  OF  THE  ROOTS 
OF  EQUATION  f{ti    )    =  0 

Let  us  examine  the  function  ^(Op)  defined  by  equation  (7^). 

Fx'oa  this  equation  we  can  deduce  that  when  9   •*•  •;-  ~,  (*  •>  -  1 

and  when  e„  *  —  or  -i-,  <2>»  +  ~.   The  sign  of  the.  derivative 
s-       k$         k 

—<-  gives  useful  information  about  the  function  fib    ).   From 

2 
equation  (7^).  we  obtain 


Af         "  9ZR  KR  kB     2  *U  %  V  ,„„ 

2   (kB  e2  -  l)    (kc  e2  -  i) 

The  abscissa  of  the  stationary  points  of  function  cKo  )  are 

given  by 

i 

ASP.      n  ram 

d^  =  °  l88J 

Employing  the  expression  given  by  equation  (87)  and  rearranging 
yields 

fkC  °2  "  V     ZU  KU  kC 
B   2  ERKRkB 


Let 


zn  K„  k*  k  0*  _  1  3 

~  -   -   r         and     ( - )  m  & 

ZR  KR  kB  kB  e   -  1 


Then 
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,zu  Ku  kc  1/3 


Solving  equation  (88)  for  stationary  points  of  function 
</>(62)  gives 

0      +   r  =  0  (89) 

which  may  be  written  as 

(<S>  +  f  )i  ®z  -  p®  +  ,» " )  =  o 

This  equation  breaks  down  into 

<Sf  ■*    P   =  °  (90) 

and 

<3>  2  -  p©  +  f  s  =  o      -  (91) 

Let  us  first  consider  equation  (86).   It  is  quadratic  with 
respect  to  Q>  .   Its  discriminant  is 

p2   _  i|f>2  .  _  3f»  <  o 

Therefore,  equation  (91)  has  no  real  roots.   Solving  equation 
(90)  yields 

<2>=  -£~i —  -  -  f 


Hence 


B2  =  82m  =  kJ-V?¥-  (92) 
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Let  us  compare  9.-   to  — . 
2m    KB 


2m  ~  kB  kc 


v 
But  _C  is  less  than  1.   Consequently 
k 

B 

K    >  i 

2H   kB    • 


Now,  6,,  can  be  written  as 
2m 


6        _  JL   (__i_t_£_) 

c    i  +  pe 

which  shows  lnimed  lately  that 


2m  <  k7 


Therefore,  we  have 


kB   e2m   kc 


It  can  thus  be  concluded  that  derivative  ~£-   vanishes  once 

0.9  o 

in  the  interval  (-  «.,  +  «,)    for  9„  =  8„  .   In  addition,  9   lies 

Z    ZJ»  2m 

between  ~  and  g~.   Equation  (53)  gives  the  following  additional 


,B 


"C 


information. 

d  <f> 


dT2  -  0 


^ne2»-L±c,   «jj.» 


8t 


when  6   •>  -—  +  0,   $J£  ->  +  » 

2   kc  -      a02 

With  this  information  we  can  set  up  the  table  for  the  variation 

of  function  f\9    )    (see  Table  1).   The  numerical  values  which 

appear  in  the  second  part  of  Table  1  are  computed  from  the 

numerical  data  given  in  Paragraph  (IV.  2).   6„  turns  out  to  be 

2m 

the  abscissa  of  a  minimum  for  the  graph  of  function  ^(8,)  (see 
Figure  19). 

e      1+  t 

2m  -  kc  +  fkB 


Thus 


r    1.50  l.ooMo.5o' 
=  ^0.375 

e  =  0.720 


,  1  +  0.720     _  1.720 

2m  B  0.25  +  0.720  x  0.50  =  OTSlO 


L  =  +  2.82 

2m 


and  the  value  of  function  <P  is  computed  as  follows  j 

f<*      )  .  1-50  x  1.00  x   0.50 
'   2M    (0.50  x  2.82  -  l)2 

1.50  x  1.50  x  0.25 

+  — * -" i" 

(0.25  x  2.82  -  1) 

--  O.lbS  +  0.087 
=  4.^6  +  6.^6  -  1 
fte2m)  „  +  9.92. 
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TABLE  1 
VARIATIONS  OF  FUNCTION  f(6  ) 


xi            en>              k- 

dej 

0.              +                +   ■> 

-     oo        -         0        +            +    «> 

-    co               -                     0 

f 

-1           """"^ 

+     CO                                             +     «• 

T  '    m' 

^^*-    -  1 

e^                                                 2                2.82                      ^                            +  «. 

AS. 

de2 

0                  +                         +    co 

-     co         -         0            +            +     oo 

-     CO             -                              0 

f 

+      CO 

-1 

+     CO                                                      +     oo 

\    / 

+9.92 

+       CO 

-1 

C6 


10. 


-0- 


I 

-2.       -0.3760  2       2.82      4 


Fig.  IS)    Function      a  ifi-  ) 
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Function  £>(9  )  is  shown  on  Figure  19.  Its  graph  intersects  the 
G  -  axis  at  two  points.  Therefore,  equation  (73)  has  two  real 
roots.   The  negative  root  is 

e2  =  -  0.376 

The  positive  root  is 

62  =  +  7.183 

Although  8  m   -0.376  is  the  ordinate  of  two  stationary  points  on 
surface  (  l)  ,  this  value  may  be  discarded  for  the  purpose  of 
optinlzation  because  it  is  negative. 

Let  us  now  check  the  two  conditions  of  feasibility, 
inequalities  (76)  and  (78).   The  first  one  is 

2S  kA  VKS^°        '  (76) 

Substituting  numerical  values  gives 

1.  x  6  x  (10. f    -   0.50  =  +  599.5  >  0 

The  second  one  is 

v      v      v  e  X      k 

"B****        ZV*V*C         „  !  >   0  (y8) 


(kB  ew3  -  i)2   (kc  ew3  -  1) 


At  the  first  stage,  6  =  0.455.   Thus,  from  equation  (12) 

si  -  o^HV.~  =  °-289 


From  equation  (17),  we  have 
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U     =   0.20(10.0   -   0.289)   =   1*9*2 
Froia   equation   (59) 

ew3  "  ofe  (1  +  rf#  ^.09 

From  equation  (7^) 

-(7.09)  o  1--5QX  1.00  x  1.50  +  1.50  x  1.50  x  0.25  _  -L 
(7.09  x  0.5  _  l)E  ,  (7.09  x  0.25  -  1)£ 

=  0.115  +  0.9^  -  1 

f(7.0Q)  =  O.059  >  0 

The  two  conditions  for  feasibility  are  satisfied  at  operating 
point  A. 
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APPENDIX  III 


NATURE  OF  THE  STATIONARY  POINTS 
OF  RESPONSE  SURFACE  (E) 


The  first  coordinate  of  stationary  points  is  given  by 
equation  (71) 


A 
Substituting  numerical  values  for  the  parameters  yields 

e  =  +  0.455 

8,  »  -  0.122 

In  addition,  equation  (73)  gives  the  following  two  values  for  9„ 

e2  =  +  7.18 

»2  «  -  0.376 

Therefore,  response  surface  (£)  has  four  stationary  points 
given  below i 

At*  0.*55.  +  7-18) 

B(-  0.122,  +  7.18) 

C(-  0.122,  -  0.376) 

D(+  O.i+55,  -  0.376) 

The  nature  of  each  of  these  stationary  points  villi  be  examined 
by  evaluation  of  the  Hessian  matrix.   As  a  general  rule,  the 
necessary  and  sufficient  condition  that  a  square  matrix  be 
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positive  definite  is  that  each  of  the  principal  minors  of  this 
matrix  be  greater  than  zero.   A  necessary  and  sufficient  condi- 
tion for  a  square  matrix  to  be  negative  definite  is  that  the 
signs  of  its  principal  minors  be  alternatively  negative  and 
positive;  that  is  the  principal  minor  of  rank  r  is  negative  if 
r  is  odd  and  positive  if  r  is  even.   We  shall  now  test  the 
matrix  H  at  the  stationary  points,  namely  A,  B,  C  and  D. 

The  general  expression  of  the  Hessian  matrix  of  the  objective 
function  J  is 


Let 


H  = 


nll 


i  a  j 

2 <*«x)" 

2  yel3e2 


1  ,!,J 

2  a^y 


**j 


2  ve,? 


1"2 


1  _^_J_ 
(262) 


S     J 


1  32J 


12  -  2  i01>S2   "  2  29"  39 


*22 


1  »    J  _ 

2T3e2)3 


From  equation  (6?),  we  have 


2_J 

■Ve " 


ZS  KS  kA 


(kA  »x  -  i) 

Dlf f erentlating  this  again  with  respect  to  9  ,  we  obtain 
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Hence 

h 


11  =  (kA  ex  -  i)= 


From  equation  ( 67 ) .  we  have 


2h  m   i  -  ZR  KR  kB zu  Ku  kc 

*e2      (kB  e2  -  i)s   (kc  e2  -  i): 


Taking  the  derivative  of  this  expression  again  with  respect  to 
eo  yields 

/*2   _   2  ZR  KR  kB  +   2  ZU  KU  *g 

^°2)2  "  (kB  62  "  1)3    (kC  62  "  1> 

and 

zR  KR  k|        Z[J  K„  k^ 


h22 


In  addition 

2, 

3.  J2 


<kB  62  "  1»3    (kC  92  "  1)3 


Henc? 


*V°2  "  °" 


h12  "  h21  -  °' 


Nov/,  the  general  Hessian  natrix  at  point  (0-,i  6„)  is 
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H  = 


(kA  6a  -  1) 


(93) 


ZR  KH  kB       ZU  KU  kC 


(kg  e2  -  l)    (kc  e2  -  i) 


Point  A  (0.^55,  7.183) 

Substituting  the  coordinates  of  point  A  in  equation  (93) 
leads  to 


H  = 
A 


1.0  x  0.50  x  (6.)a 
(6.00  x  0.^55  -  l)3 


1.  50x1. 00x( 0.  50 ?3  1.  50x1  .50x(0.?-5) 
(0.50x7.183  -  l)3  (0.25x7.183  -  X) 


HA  = 


18 


(1.73) 


Finally 


HA  = 


+3.W 
o 


0 

_Oi27JL 

♦ 

0 

UK) 

(2.59)3 

(0 

796)3 

i, 

+  0.299 


H  is  positive  definite.   The  two  principal  minors  are  +3.^8  >  0 
A 

and  (+  2.48X+  0.299)  -  0  >  0.   Hence  point  A(+  0.^55,  +  7.183) 
is  a  minimum  en  response  surface  (£).   The  procedure  set  up  in 
Section  IV. 5  to  generate  the  contours  around  this  point  has  been 
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implemented  by.  means  of  subroutine  ISOBA'T  (see  Figure  20). 

Point  B  (-  0.122,  +  7.18) 

Substituting  e^  =  -  0.122  end.  62  ■  +  7.18  in  to  equation  (88) 


yields 


18 


t -1.732 r 

0 

-  3-46 
0 


,  Q.J15__  .     o.i^o 
(2.59)3   (0.796)3 


0 
+  0.299 


H  is  not  positive  nor  negative  definite  because  Its  principal 

s 

minors  are  -  3.46  <  0  and  (-  3.46) (+  0.299)  <  0.  As  a  result, 
point  B(-  0.122,  +  7.18)  is  a  saddle  point  on  response  surface 
fc). 

Point  C  (-  0.122,-  0.376) 

This  point  has   the  abscissa  of  B  and  the   ordinate   9     = 

2 

-   0.3?6.      Thus 


hll  "    "   ?•* 


and 


22 


0.375 


0.140 


[0.50   x   (-0.376)    -  1]  [0.25   x   (-0.376)    -  i; 


0.375  O.lftO 

(-1.188)  (-1.094)" 


9k 


<>**        ) 

Specifications 

3=( 

AINT(10.*S;  r;)+i.)/ir>. 

3110=0. 

Do 

110     1=1,10 

y={m.  r(r;:irug),.io.)-;-i.)/ij.  | 

Zo   30      J=l,25 


\j.j?;      ti,x2,::  / 


^L 


"©    .        UEjHB 


EM> 


-'© 


Figure  20(a).  Flow  -  chart  of  subroutine  ISOBAT 
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T 


|    Y=(AIHT(Ti-:iK(2)>10.))AO-. 

Do  po     Jsl,20 


<Y-l/lf3   <10" 
or 
5C-1/KC   < 


<C 


3110   =    1. 

3=i.io..s;:i 


[S 


Ye-- 


Yes 


-CMID 


Figure     20(b).  Bad  of  flow  -  chart  of  subroutine  TSORAT 
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=  -  0.331 
The  Hessian  matrix  at  point  C  Is 
-  3.46        0 
0       -  0.331 


:V 


H  is  negative  definite  because  its  principal  minors  are 

-  3.46  <  0  and  (-3.46) (-0.331)  >  0.   Therefore,  point  C  (-0.122, 

-0.376)  is  a  maximum  on  (£). 

Point  D  (+0.455,-0.376) 

This  point  has  the  same  abscissa  as  A  and  the  same  ordinate 
as  C.   Matrix  E  can  immediately  be  written  as 


+  3.48 


0 
-0.331 


H  is  neither  positive  definite  nor  negative  definite  because 
its  principal  minors  are  +3.48  >  0  and  (+3.48) (-0. 331)  <  0. 
Therefore,  point  D  is  a  saddle  point  on  response  surface  (E). 
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APPENDIX  IV 

COMPUTATIONAL  SCHEME  FOR  THE 
SEARCH  TECHNIQUE  APPROACH 

Basic  features  of  the  sequential  Simplex  technique  (11) • 

We  consider,  initially,  the  minimization  of  a  function  of  n 

variables,  without  constraints.   P.,  P, ,  ...,  P  are  the  (n+1) 

points  of  the  current  simplex  in  n-dimensiohal  space.   We  write 

y  for  the  function  value  at  P  and  define  h  as  the  suffix  such 

that 

y  =  ma-(y  ) 
h    t   i 

and  £  as  the  suffix  such  that 

y-  =  min  (y  ) 
*     i    1 

Further  we  define  P  as  the  centroid  of  the  points  with  i  f.   h, 

and  write  [P,  P.,1  for  the  distance  from  P,  to  P..   At  each  stage 
L  i   j  i     j  e 

in  the  process  P  is  replaced  by  a  new  point;  three  operations 
are  used  -  REFLECTION,  CONTRACTION,  and  EXPANSION.   These  are 
defined  as  follows i   The  reflection  of  P  is  denoted  by  P*f  and 
its  co-ordinates  are  defined  by  the  relation 

P*  =  (1  +  ol)P  -  aPh 

where  a  Is  a  positive  constant,  the  reflection  coefficient. 
Thus  F*  is  on  the  line  joining  F  and  F,  on  the  far  side  of  P 
from  Ph  with  [p*  f]  .  a[Ph  P],   If  y*   lies  between  y  and  y„, 
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then 

Ph 

is  raplac 

:ed  by  F*  and  we  start  again  with  ' 

the  new  simplex. 

If 

y  <  y(,  i 

. .e.,  if  reflection  has  produc 

:ed  a 

new  minimum, 

then 

we 

expand  P* 

■  to  P»*  by  the  relation 

P* <;  a  YP* 

+  (i  -  y)p 

The 

expansion  coefficient  y>  which  is  greater 

than 

unity,  is 

the 

rat; 

i.o  of  the 

distance  [  P**  P]  to  [  P*  f]  . 

If  y 

»«  j.   y  We 

repli 

ace 

P  by  P«* 
h 

■   and  restart  the  process ;   but 

!  if  , 

F**  >  yp  then 

we  have 

failed  to  find  a  better  point  by  expansion 

,  and  we 

repli 

ace 

Ph  by  F* 

before  restarting. 

If  on  reflecting  P  to  P*  we  find  that  y*  >  y   for  all 

i  =  h,  i.e.  that  replacing  P,  by  F*  leaves  y*  the  maximum,  then 
/  h 

we  define  a  new  P  to  be  either  the  old  P,  or  P* ,  whichever  has 
h  h 

the  lower  y  value,  and  form 

P*«  =  8P  +  (1  -  p)P. 

The  contraction  coefficient  0  lies  between  0  and  1  and  la  the 

ratio  of  the  distance  [P**  P]  to  [P  T],    We  then  accept  ?**  for 

P  and  restart,  unless  y9*  >  min(y  ,  y*),  i.e.,  the  contracted 

point  is  worse  than  the  better  of  F  and  P*.   For  such  a  failed 

h 

contraction  we  replace  all  the  P.'s  by  (P  +  ?_^)/2  and  restart 
the  process. 

The  criterion  adopted  to  stop  the  search  is  to  compare  the 


"standard  error"  of  the  y's  in  the  form  jz{y     -  y)s/n  with  a 
pre-set  value. 

The  bcgm  SIMPLE  has  been  given  to  the  corresponding  FOflTHAM 
subroutine  for  this  search  procedure. 
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In  order  to  perform  a  search  the  value  of  the  objective 
function  should  be  available  for  a  given  decision  vector.  The 
FORTRAN  name  of  this  subroutine  is  OBJ.   Flow-chart  of  subroutine 
SIMPLE  appears  in  Figure  21. 

Constraints  on  the  volume  to  be  searched . 

If,  for  example,  one  of  the  x  must  be  non-negative  in  a 
minimization  problem,  then  the  original  sequential  simplex 
search  method  may  be  adopted  as  follows!   The  scale  of  the  x 
concerned  can  be  transformed,  e.g.,  the  function  can  be  modified 
to  take  a  large  positive  value  for  all  negative  x .   In  the  latter 
case  any  trespassing  by  the  simplex  over  the  border  villi  be 
followed  automatically  by  contraction  moves  which  will  eventually 
keep  it  Inside.  TM  r,   ftsthod  is  Illustrated  in  Figure  22.   If 
the  reflected  point  P*  trespasses  the  border  line  of  the  permitted 
region,  a  contraction  follows,  which  results  in  the  new  simple:': 
V   ,  Pp,  F**.   If  the  expanded  point  Pft*  trespasses  the  border 
line,  P**  is  replaced  by  the  reflected  point  P»  which  was  within 
the  permitted  region.   In  either  case  the  new  simplex  is  within 
the  constraint  and  the  iteration  process  can  be  carried  on. 

Determination  of  the  optimal  policy  by  the  empirical  search . 
The  various  components  of  the  optimization  problem  have 
been  stated  in  Section  IV.   The  procedure  consists  of  the  follow- 
ing steps t 

Compute  9,    ,  by  equation  (72)  and  check  the  first 
1  opt 

condition,  equation  (76). 
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Specification 


.  IT  SH  =   0 
ICONV  =   0 


Initial  functii 
values 


)eflne  the  woi 
point    r<JH) 


Define   the   oe 
point  T(JL) 


Figure  21  fa).     Plow  -  chart  of  subroutine  SBTPUE 
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Yet 


For  definition 
of  eentroid ,    all  vertices 
are  equally  weighted 


The  best  point 
has   the   largest 
weight 


Coaput e   coord  i na t  e s 
of  REFLECTED  POI1IT   : 


L 150        COl'TI    U5 


Yes 


Figure  2Kb).    Continuation 
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Conpute  coordinates 

of  CONTHACTSC  POIOT  SCO]: 


I 


Figure  21(c).    Continuation 
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(E> 


Yes 


Replace  the  worst 
point  by  contracts 
point 


Contract   the  whole 
Simpler   towards   the 
best   point 


Joapute  C   I  1 
JITE:w:TIT2K+l 


Figure  21(d) .     End  of  subroutine  SIMPLE 


10't 


IlQs^_oJL_refLGct:on 


—  Old   simplex 
New  simplex 


Cass 


Qi     expansion 


Fig. 22.  Treatment    of    inequality    consents    on    the    ind°D«nd*n 


vanao!;j3     in     simplex    technique 
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Determine  9_   .  by  solving  equation  (?3).   Compute  0  .  and 
£■    opt  w2 

check  the  second  condition,  equation  (?8). 
Step  3 

Initialize  counters  and  set  the  first  guess  for  optimal 
solution. 
Step  h 

In  order  to  perform  the  search  an  Initial  simplex  and  the 
inequality  constraints  should  be  specified.   Thus  this  step 
consists  of  sotting  the  initial  simplex  in  the  neighborhood  of 
the  optimal  point  already  defined  in  Steps  1  and  2  and  computing 
corresponding  to  the  initial  guess  for  9„  set  in  Step  J. 


Call  the  search  technique.   If  the  resulting  point  is  close 
enough  to  the  initial  guess,  its  coordinates  constitute  the 
optimal  policy.   If  not,  the  result  of  the  search  serves  as 
starting  point  for  the  next  iteration. 
Step__6 

After  convergence  the  corresponding  state  variables  can  be 
computed  using  the  performance  equations. 
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APPENDIX  V 
DESCRIPTION  OF  MAIN  PROGRAM  MICU20 

The  main  program  for  solving  this  problem  has  been  given 
the  name  MICU20.   In  this  paragraph  all  FORTRAN  locations  will  be 
written  with  capital  characters. 

The  program  MICU20  consists  of  six  parts  as  illustrated  in 
Figure  ?3. 
Fi_r_st_ga r  1 1 

It  contains  the  specif fications,  the  definition  of  statement 
functions  needed  in  the  remainder  of  the  program,  the  input 
statements  for  the  physical  and  economical  parameters  and  the 
input  statements  of  option  parameters  NOPT,  ICURV  and  ICONT  and 

the  search  parameters.   In  addition,  it  performs  the  echo-check 

i 
of  these  data. 

Second  parti 

It  is  devoted  to  the  determination  of  the  optimal  policy. 
Note  that  equation  (73)  which  is  of  degree  4,  has  been  solved  by 
the  NEWTON  method.   The  name  of  the  subroutine  which  performs 
this  search  is  precisely  NEWTON. 

The  values  of  the  first  members  of  inequalities  given  by 
equations  (76)  and  (78)  are  called  respectively  TEST1  and  TEST2. 
Th_irfi_parti  •'  ' 

This  part  deals  with  the  case  where  the  operating  policy  (9-j, 
8?)  Is  fixed  by  user  without  being  optimal. 

It  simply  consists  of  reading  the  chosen  decision  variables, 
and  calling  subroutine  OBJ  which  gives  the  value  of  the  objective 
function. 


CUnz) 


10? 


Specifications 


[  Define  statement   functions 

\  BEAD/.raiTE  7 

\  XISO.KA.KS  ,  ZAS  ,  JTRS.rUS  / 
\     K3,KR,YB  t.YPH  / 

\   :;c,ku,xcu,::pu  / 


READ/  ./..I  ?E 

ALPHA, 3SIA,  fA!  ..A 

STEP.IHAJBX 


HEAD/     IT 
gpsTL   r  T.  I 


Figure  23(a).     Flow-chart  of  cain  program  UICU20 

First  part:  Specifications,   statenient  functions 
&  data  echo  cheek 
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X 


Co  route  Ti  10?T,    rESTl 
CALL  NBWTOH    (TH20PT) 
Compute  r  .  .'3 
TSST2,    QUA.HT(TH'J3) 


\= 


XE10PT,  TE3T1,  IH20PT, 


Initialize  coun 
Sot  Initial  gue 


S3     ,"n   1.,     PH2 


\    Set  Initial  simplex  &  inequa-  / 
\ 11  ty  constraints  / 


TWIN(l)-  LTH1  <ETA 

X  'J 
:■:  I-  (2)-KTH2   <S?A 


G3D 


Figure  23(b).     Second  part:     Determination  of  optimal  oolicy 
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il  =   E-!IW(1) 


EDI  =  I«H{1) 
KT2  =   T:_I/(2) 
Confute  r6;».T 


\  WRITS  LL.iC  , !  :?1 ,  HT2  ,i.rl 


Figure  23(c).     End  of  second  part 
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\±Ei;ii±±/ 


DV(1)   =   HT1 

DV(2)   =   BT2 

Compute  TH72'i  EHW3 


CALL   0    ,?(DV,    3) 


Figure  23(d).  Third  part:  Case  of  non-cptinal  policy 


Hi 


\'m"is  ^-^7 


ail  >  ra 

-\ 

b 

t  first 

variab: 
stage 

-/ 

\5 

ITS  TE.^2 

,Ta  0. 

-/ 

Arcm-m*./ 


\  '.;,;IJ.ri  sts'.Se  v-ri -:-]"'.  r:  7 
\         at   second   stage        / 


Figure     23(e).     Fourth  part:     State  variables   of  optinal  policy 
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Figure  23(f).     Fifth  part:     Concentration  curves 
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\hbad  delts,  deltax/ 


CALL   IS03AT 


•KOHIiAL  TEl&aKATIOtJ 


Figure  23(g).     Part  5]:     Drawing  of  contours 
End  of  program  ?TICU20 


11'+ 


Z?il£t h  part  i 

The  state  variables  resulting  from  the  previously  chosen 
policy  are  determined. 
Fifth  cart  I 

In  this  part  the  concentration  curves  are  plotted.   In 
addition  the  slopes  at  the  points  of  discontinuity  are  determined, 
according  to  Appendix  I. 
Sixth  parti 

By  calling  subroutine  ISOBAT  various  contours  around  the 
optimal  point  are  generated  using  the  procedure  described  in 
section  IV. 5.   The  meaning  of  the  symbols  in  Figure  23  appears 
in  Table  2.   Table  3  contains  the  list  of  the  required  data  cards. 
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TABLE  2 
MAIN  FORTRAN  SYMBOLS 


Symbol 


Meaning  or  corresponding  algebraic  variable 


MICU20 
SIMPLE 
NEWTON 

ISOBAT 

OBJ 
QUART 

KA,  KB,  KC 
KS,  KR,  KU 

ma,  m'jr,  kuu 

THW1,  THW2,  THW3 
DELTAS,  DELTAI 

NDIM 

ALPHA,  BETA,  GAHHi 

IBARI 

EPSIL 

LIMIT 

ICONV 

TMIM 
NITER 


Main  routine 

Subroutine  for  SIMPLEX  search  technique 

Subroutine  implementing  NEWTON'S  method  for 
root  searching 

Subroutine  to  generate  the  contours  around 
the  optimal  point  on  the  response  surface 

Subroutine  for  computing  objective  function 

Function  f(Q?) 


Kg,  KH,  KTJ 


Hwl'  HW2'  Uw3 

Increment  on  objective  function  and  0„ 

respectively  for  contours  drawing 

Dimensionality  of  search 

a,  3,  Y 

Option  parameter  for  SIMPLE 

Convergence  criterion  for  SIMPLE 

Maximum  number  of  iterations  to  be  performed 
during  the  search 

Convergence  index 

Optimal  decision  vector 

Number  of  iterations  actually  performed  in 
SIMPLE 


HEVAI 


Number  of  function  evaluations 
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TABLE  2  (Cont.) 


X1S1,  X2A1,  X1R1, 
X1U1 

DELTAS 

DELTAR 

DELTAU 

X1R2,  X1U2,  X2B2, 
• X2C2,  X1P2 

R,  BET 

W2,  W3 

NOPT,    1CURV, 
ICONS 

ETA 

MAXIT 

TEST1 

TE3T2 

TK10PT,  TH20PT' 

LLAC 

TH 

ET1,  RT2 

S 

DX131 

DX2A1 
DX1R1,  DX1U1 


Si(  A,  R1,  0, 


S0  "  sl 


R2.  U2, 


i,    C,  P 


Recycle  ratio  and  clarlfyer  efficiency 

ew2'  ew3 

Option  parameters 

Criterion  for  over-all  convergence 

Maximum  number  of  iterations  for  optimisation 
loop 

Value  of  first  member  of  inequality  (76) 

Value  of  first  member  of  inequality  (?8) 

eiopt'  62opt 
Iterations  counter 

Current  Simplex 

Operating  residence  times 

Objective  function 

dei  I  e*i 

4JLI 

3i|e^ 


de, 


dRj^ 

aw: 


dU, 


TABLE  2  (Cont.) 
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T 

SH 

SI, 

TBAH 

SBAR 

TREF, 

SREF 

TEXP , 

SEXP 

Current  decision  vector  in  SIKPLE 

Highest  function  value  (y,  ) 
h 

Lowest  function  value  (y .) 

Coordinates  of  the  centroid 

Value  of  the  objective  function  at  centroid 

Coordinates  and  function  value  at  reflected 
point 

Coordinates  and  function  value  at  expanded 
point 
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Card 


TABLE  3 
DATA  CARDS 


I  and  2  Title 

3  JOSO i  EA,  KS,  YAS,  YRS,  YUS 

k  KB,  KB,  YBR,  YPR 

5  KC,  KU,  YCU,  YPU 

6  MUS,  KUR,  MUU 

7  NOPT,  ICURV,  ICONT 

8  R,  BET 

*  If  NOPT  =  0  and  ICONT  =  0 

9  RTl,  RT2 

*  If  NOPT  ^  0  and  ICOMT  =  0 

10  ALPHA,  BETA,  GAMMA,  STEP,  IBARY 

II  EPSIL,  LIMIT,  ETA,  KAXIT 

*  If  NOPT  /  0  and  ICONT  ^  0 

12  DELTS,  DELTAY 
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ABSTRACT 

This  work  considers  an  anaerobic  digestion  process  represented 
by  a  two-step  mechanism  and  carried  out  in  a  two-stage  continuous 
digester  system.   A  mathematical  formulation  cf  the  kinetic  and 
flow  models  has  been  given,  which  allows  simulation  of  the  system. 
On  this  "oasis,  an  objective  function  of  economic  type  can  be 
constructed,  the  minimization  of  which  yields  the  optimal  design. 
The  analysis  of  the  system  by  means  of  this  mathematical 
formulation  also  shows  that  wash-out  times  are  important  design 
parameters.   Finally,  a  contact  process  in  which  organisms  are 
recycled  to  the  second  stage  is  compared  to  the  conventional 
process.   From  this  comparison  it  can  be  concluded  that  the 
contact  process  has  definite  advantages  over  the  conventional 
process.  t 


